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Abstract 

We review the Resonating Valence Bond (RVB) theory of high temper- 
ature superconductivity using Gutzwiller projected wave functions that 
incorporate strong correlations. After a general overview of the phe- 
nomenon of high temperature superconductivity, we discuss Anderson's 
RVB picture and its implementation by renormalised mean field theory 
(RMFT) and variational Monte Carlo (VMC) techniques. We review 
RMFT and VMC results with an emphasis on recent developments in 
extending VMC and RMFT techniques to excited states. We compare re- 
sults obtained from these methods with angle resolved photoemission spec- 
troscopy (ARPES) and scanning tunnelling microscopy (STM). We con- 
clude by summarising recent successes of this approach and discuss open 
problems that need to be solved for a consistent and complete description 
of high temperature superconductivity using Gutzwiller projected wave 
functions. 
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1 Introduction 

This paper reviews developments in the use of Gutzwillcr projected wave func- 
tions and the Resonating Valence Bond (RVB) theory in the context of high 
temperature superconductivity. We attempt to review comprehensively, both 
the general framework of the Gutzwiller-RVB theory and to summarise several 
recent results in this field. Though many of these results were indeed motivated 
by the phenomenon of high temperature superconductivity and the rich phase 
diagram of these compounds, it is not our intention to review high temperature 
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Figure 1: Crystal structure of La2Cu04. Left panel shows the layer structure 
along the c-axis, the right panel the structure of the Cu02 plane. From [2]. 



superconductivity per se. Nonetheless, it is well nigh impossible, if not meaning- 
less, to attempt to write a review of this nature without discussing certain key 
experimental results. Our choice in this matter is dictated by the fact that most 
techniques used in the study of Gutzwiller projected wave functions address the 
calculation of single particle spectral features. Consequently, after discussing 
some basic facts and a historical perspective of the Gutzwiller- RVB concept, we 
present an overview of experimental results from angle resolved photoemission 
spectroscopy (ARPES) and scanning tunnelling microscopy (STM) within this 
introductory section. We also discuss briefly, a few alternative theories based 
on repulsive electronic models, to illustrate the complexity of the subject. 

1.1 High temperature superconductivity 

Twenty years ago Bednorz and Miiller [T] discovered high temperature super- 
conductivity in Sr-doped La2Cu04. Subsequently high temperature supercon- 
ductivity was reported in many other Cuprates. These compounds have a lay- 
ered structure made up of one or more copper-oxygen planes (see figure [T]) . It 
was soon realised that many of the HTSC have an insulating antiferromagnetic 
parent compound that becomes superconducting when doped with holes or elec- 
trons. This is fundamentally different from, say, superconductivity in alkaline 
metals and clearly calls for a novel mechanism. 

These unusual observations stimulated an enormous amount of experimen- 
tal as well as theoretical works on HTSC, which brought about numerous new 
insights into these fascinating compounds. The d-wave nature of the super- 
conducting pairs [3] as well as the generic temperature-doping phase diagram 
(figure[2]) are now well established. On the theoretical front, several approaches 
successfully describe at least some features of the HTSC. In addition, new sophis- 
ticated numerical techniques provide us with a better understanding of strong 
correlation effects that are clearly present in the HTSC. Progress in the field 
of high temperature superconductivity has also influenced many other flelds in 
condensed matter physics greatly. Research on HTSC has a very fruitful history, 
and continues to broaden our knowledge of strongly correlated electron systems. 
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Given the numerous theories advanced to explain the phenomenon of high 
temperature superconductivity [4], it is important to examine carefully the 
strengths and weaknesses of any given theoretical approach and its relevance 
to experimental observations. In this review, we examine the resonating valence 
bond (RVB) scenario which proposes a simple, yet nontrivial wave function to 
describe the ground state of Mott Hubbard superconductors, i.e. superconduc- 
tors that are obtained by doping a Mott Hubbard insulator. Wc discuss various 
theoretical calculations based on the so called Gutzwiller-RVB wave function 
both in the context of our work [SJ El [HI [SI HH] and other recent developments. 

The Gutzwiller-RVB theory provides a direct description of strongly corre- 
lated superconductors. An advantage of this approach is that the theory can be 
studied by a variety of approximate analytical techniques as well as numerical 
methods. We will discuss later how the theory yields many results that are 
in broad agreement with various key experimental facts. However, to obtain a 
more complete description of HTSC, the Gutzwiller-RVB calculations need to be 
be extended to be able to describe finite temperature and dynamic effects. This 
review should provide an adequate starting point for further extensions of this 
method as well as phenomenological calculations of various physical quantities 
that are relevant to the phenomenon of high temperature superconductivity. 

1.2 A historical perspective 

The notion of Resonating Valence Bonds was introduced by Pauling [TTl [T^] 
in the context of the Heitlcr-London approximation for certain types of non- 
classical molecular structures. Anderson and Fazckas [TSl [IH then generalised 
this concept to the case of frustrated magnetism of localised spin- 1/2 mo- 
ments. The RVB theory came to a first full bloom with the discovery of high- 
temperature superconductivity when Anderson [15j suggested that an RVB state 
naturally leads to incipient superconductivity from preformed singlet pairs in 
the parent insulating state. 

A detailed account of the progress made after Anderson's seminal RVB pro- 
posal will be presented in this review in subsequent sections. At this point we 
will make a few comments regarding the general lines of development of the 
theory. 

The core of the RVB concept is variational in nature; the RVB state may be 
regarded as an unstable fixed point leading to various instabilities, such as an- 
tiferromagnetic order, superconductivity, etc., very much like the Fermi-liquid 
state. However, in contrast to Fermi-liquid theory, there is no simple Hamilto- 
nian known for which the RVB states discussed in this review are exact solutions. 
For this reason, the theory developed historically along several complementary 
lines. The first one is the quantification of the variational approach via the 
variational Monte Carlo method (VMC). This approach was initially hampered 
by the problem of implementing the numerical evaluation of a general RVB 
wave function algorithmically |16j . But when this problem was solved |17j , the 
method evolved quickly into a standard numerical technique. 

Very early on it was realized [TH] , that essential aspects of the RVB concept 
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could be formulated within a slave-boson approach, which led to the develop- 
ment of gauge theories for strongly correlated electronic systems in general, and 
high temperature superconductivity in particular. This line of thought has been 
reviewed comprehensively by Lee, Nagaosa, and Wen [T^ . 

The superconducting state is an ordered state and this statement applies also 
to the case of the high temperature superconductors. Where there is an order 
parameter, there is a mean field and it was felt early on that a suitable mean- field 
theory should be possible when formulated in the correct Hilbert space, using 
the appropriate order parameters. This line of thought led to the development 
of the renormalised mean-field theory (RMFT) [5D]. This theory will play a 
prominent role in this review, as it allows for qualitative analytical predictions 
and, in some cases, also for quantitative evaluations of experimentally accessible 
response functions. 

There is a certain historical oddity concerning the development of the RVB 
concept and of the theory. After an initial flurry, there was relatively little 
activity in the 1990's and the Gutzwillcr-RVB approach returned into the cen- 
ter of scientific interest only in the last decade with the evaluation of several 
new response functions [21] , allowing for a detailed comparison with the (then) 
newly available experimental results. In retrospect, is not quite clear why this 
particular approach lay idle for nearly a decade. It is tempting to speculate that 
perhaps the concept was too successful initially, predicting d-wave superconduc- 
tivity in the cuprates at a time when available experimental results favoured an 
s-wave. 

1.3 Experiments 

The discovery of high temperature superconductivity stimulated the develop- 
ment of several new experimental techniques. Here, we shall mention some key 
experimental facts concerning the HTSC and refer the reader to more detailed 
summaries of experimental results, available in the literature p,[51 [Tni[^[^[^ . 

An early and significant result was the realization that HTSC are doped 
Mott insulators, as shown in the generic temperature-doping phase diagram 
(see figuref2|). The figure shows the antiferromagnetic phase in the undoped 



(half-fille43) CO mpound with a Nccl temperature of about T/y ~ 300K . Upon 
doping, antiferromagnetism is suppressed and superconductivity emerges. The 
behaviour of Tc with doping exhibits a characteristic "dome" . While clcctron- 
and hole-doped HTSC share many common features, they do exhibit some sig- 
nificant differences, e.g. the antiferromagnetic region persists to much higher 
doping levels for electron-doped Cuprates. 

We will restrict our attention to the hole-doped compounds, partly because 
they are better characterised and more extensively investigated, and also be- 
cause the hole-doped HTSC exhibit a so-called pseudogap phase (with a par- 
tially gapped excitation spectrum) above the superconducting dome. The onset 

^The copper ion is in a cfi configuration, with a single hole in the d-shell per unit coll. As 
shown by Zhang and Rice |25j this situation corresponds to a half-filled band in an effective 
single-band model. 
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Figure 2: Generic phase diagram for the high temperature superconductors 
(antiferromagnetic region AF, superconducting phase SC). The temperature 
below which superconductivity (a pseudogap) is observed is denoted by Tc (T*). 
T* is possibly a crossover temperature, though some experiments (compare 
figure [7]) indicate a relation to a mean-field like second order transition. 
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Figure 3: A schematic picture of the 2D Fermi surface (thick black line) of 
HTSC in the first quadrant of the first Brillouin zone. The lattice constant a is 
set to unity. The (p defines the Fermi surface angle. 
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Figure 4: Energy distribution curve (EDC) at fixed momentum k = (tt, 0) for an 
overdoped (87K) Bi2Sr2CaCu208+5 (Bi2212) sample in tlie normal state (NS) 
and superconducting state (SC). From [2]. 

temperature of the pseudogap decreases linearly with doping and disappears in 
the overdope41 regime. The origin of the pseudogap is one of the most contro- 
versial topics in the high-Tc debate. The relationship between the pseudogap 
and other important features such as the presence of a Nernst phase [26l [27] , 
charge inhomogcneities |28j . the neutron scattering resonance |29| . marginal 
Fermi liquid behaviour [30], or disorder [31]. For a detailed discussion of the 
pseudogap problem, we refer to a recent article by Norman, Pines, and Kallin 



We now discuss some results from angle resolved photoemission spectroscopy 
(ARPES) and scanning tunnelling microscopy (STM). since they are immedi- 
ately relevant to the theoretical considerations and results presented in the later 
sections. These two techniques have seen significant advances in recent years 
and provided us new insights on the nature of the pseudogap, superconducting 
gap and quasiparticles in the superconducting state. As we will show in the 
following sections, many features reported by these experiments can be well 
understood within the framework of the Gutzwiller-RVB theory. 

1.3.1 Angle resolved photoemission spectroscopy (ARPES) 

By measuring the energy and momentum of photo-electrons, ARPES provides 
information about the single particle spectral function, A(k,a;). The latter 
quantity is related to the electron Green's function by A(k,uj) = — ilmG(k, w) 
[35] . In this subsection, we summarise some key results from ARPES that any 

^The superconducting phase is often divided into an optimal doped (doping level with 
highest Tc), an overdoped (doping level higher than optimal doped), and an underdoped 
(doping level lower than optimal doped) regime. 
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Figure 5: Momentum dependence of the spectral gap A at the FS in the super- 
conducting state of an overdoped Bi2212 sample from ARPES. The black line 
is a fit to the data. For a definition of the FS angle 4> see figure [3] From [33] . 

theory of HTSC has to address. The reader is referred to the extensive ARPES 
reviews by Damascelli, et al. j23| and Campuzano, et al. |24| for a discussion on 
experimental detail. 

In figure [3] we illustrate a schematic picture of the two-dimensional (2D) 
Fermi surface (FS) of HTSC in the first quadrant of the first Brillouin zone. It 
can be obtained by ARPES scans along different angles (/>. The FS for each (j> 
is then determined in general (but not in the underdoped region [9]) by looking 
at the minimum energy of the photoelectron along this direction in momentum 
space. A typical energy distribution curve (EDC), i.e. photoemission intensity 
as a function of energy at fixed momentum, from an ARPES experiment is shown 
in figured) The figure shows the photoemission intensity at the (tt, 0)-point of 
a photoelectron in the superconducting state (T ^ T^) and in the normal state 
(Tj, > T). In the superconducting state, one sees the characteristic peak-dip- 
hump structure; the peak can be associated with a coherent quasiparticle. Above 
Tj,, coherence is lost and the sharp peak disappears. 

In the early years following the discovery of high temperature superconduc- 
tivity, it was unclear if the pairing symmetry were isotropic (s-wave like), as 
in conventional phonon-mediated superconductors, or anisotropic. Later exper- 
iments have consistently confirmed an anisotropic gap with c?-wave symmetry 
[3]. The angular dependence of the gap function is nicely seen in ARPES mea- 
surements on HTSC (figureini), which accurately determine the superconducting 
gap |Ak| at the FS. As illustrated in figure [5] for a Bi2Sr2CaCu208+5 (Bi2212) 
sample, the gap vanishes for cj) = 45°. This direction is often referred to as 
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Figure 6: Doping dependence of T* (the onset of the pseudogap, compare with 
figure [2]), and of the peak and hump binding energies in the superconducting 
state (see figure |4]). The empirical relation between Tc and doping x is given by 
TJT^^^ = 1 - 82.6(a; - 0.16)^ with T^'^'^Qb K. Data for Bi2212, from [34]. 
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Figure 7: T* /Tc versus 2A/{kBTc) for various Cuprates compared to the mean 
field relation, 2A/ (kBT*) = 4.3, valid for d-wave superconductivity [35], where 
T* replaces Tc. From [36]. 



the 'nodal direction', the point at the FS is then called the 'nodal point' or 
'Fermi point'. In contrary, the gap becomes maximal for = 0°, 90°, i.e. at 
the 'anti-nodal point'. 

Another feature well established by ARPES is the doping dependence of 
the superconducting gap and the opening of the pseudogap at a temperature 
T* > Tc. Unlike conventional superconductors, HTSC exhibit a strong deviation 
from the BCS-raticH of 2A/(fcBTc) « 4.3 for superconductors with a d-wave gap 
function. In HTSC, this ratio is strongly doping dependent and becomes quite 
large for underdoped samples, where the transition temperature Tc decreases, 
while the magnitude of the superconducting gap increases. As illustrated in 
figure [6] for a Bi2212 sample, the binding energy of the peak at (tt, 0), i.e. the 
superconducting gapQ, increases linearly (with doping) while approaching the 
half filled limit. Interestingly, the opening of the pseudogap at temperature T* 
seems to be related to the magnitude of the gap. The modified ratio 2 A/ (kBT*) 
is a constant for HTSC at all doping levels and the constant is in agreement with 
the BCS ratio, 4.3 (sec figure [7|), with Tc substituted by T* . This experimental 
result is as a remarkable confirmation of early predictions from Gutzwiller-RVB 
theory, as we will discuss in further detail in latter sections. Figure[6]also reveals 
that the hump feature (see EDC in figure [4]) scales with the binding energy of 
the peak at (tt, 0). 

An additional doping dependent feature extracted from ARPES data is the 
spectral weight of the coherent quasiparticle (QP) peak. Feng, et al. [37] defined 

■^The weak coupling BCS-ratio for s-wavc superconductors, 2A/(kgTc) ~ 3.5. 
''When speaking about (the magnitude of) the superconducting gap A in a d-wave state 
without specifying the momentum k, we mean the size of the gap | Ak| at k = (vr, 0). 
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Figure 8: a) Doping dependence of the supereonducting state spectra in Bi2212 
at (tTjO) taken at T <C Tc. The doping level is decreasing form the top curve 
downwards. Samples are denoted by OD (overdoped), OP (optimal doped), and 
UD (underdoped), respectively, together with their Tc in Kelvin, e.g. OD75 
denotes an overdoped sample with Tc ~ 75K. b) The doping dependence of 
supereonducting peak ratio (spectral weight of coherent peak with respect to 
the total spectral weight) is plotted over a typical Bi2212 phase diagram for 
the spectra in a). AF, antiferromagnetic regime; SC, superconducting regime. 
From [37]. 

Figure 9: Electron dynamics in the La2-.TSri:Cu04 (LSCO) system, a) Disper- 
sion energy, E, as a function of momentum, k, of LSCO samples with various 
dopings measured along the nodal direction. The arrow indicates the posi- 
tion of the kink that separates the dispersion into high-energy and low-energy 
parts with different slopes. Ej? and \ip, are Fermi energy and Fermi momen- 
tum, respectively, b) Scattering rate as measured by MDC width of the LSCO 
(x=0.063). From [39] . 

a superconducting peak ration (SPR) by comparing the area under the coherent 
peak with that of the total spectral weight. Figure [5] depicts EDCs at several 
doping levels together with the computed SPR as a function of doping. The QP 
spectral weight strongly decreases with decreasing doping and finally vanishes 
[37l[38]. Such a behaviour is well understood by invoking the projected nature 
of the superconducting state as we will discuss in the following sections. 

Since ARPES is both a momentum and energy resolved probe, it allows for 
the measurement of the dispersion of the coherent peak. Here, we concentrate on 
the nodal point, where the excitations are gapless even in the superconducting 
state, owing to the d-wave symmetry of the gap. The dispersion around the 
nodal point is well approximated by Dirac cones, whose shape is characterised 
by two velocity, vp and v^. The Fermi velocity vp is determined by the slope 
of the dispersion along the nodal direction at the nodal point, whereas the gap 
velocity va is defined by the slope of the 'dispersion' perpendicular to the nodal 
direction at the nodal point. Since all other k-points are gapped, the shape of 
the Dirac-like dispersion around the nodal point is of particular importance for 
the description of any effect depending on low- lying excitations. 

FigurclHKa) illustrates the slope of the dispersion along the nodal direction for 
La2-a;Sra;Cu04 (LSCO) samples at various dopings. The ARPES data reveals 
a significant splitting in high-energy and low-energy parts, whereas the low- 
energy part corresponds to the Fermi velocity vp. Within ARPES data [see 
figure IHlJa)] the Fermi velocity vp is only weakly doping-dependent. ARPES 
can also determine the gap velocity va by looking at the spectral gap along the 
Fermi surface as done in figure [5] Together with the vp, the va determines 
the shape of the Dirac cones, which, according to ARPES, is quite anisotropic 
(vp/vA ~ 20 around optimal doping) [33j . This result is confirmed by thermal 
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Figure 10: Schematic illustration of the temperature evolution of the Fermi 
surface in underdoped Cuprates as observed by ARPES. The d-wave node below 
Tc (left panel) becomes a gapless arc above Tc (middle panel) which expands 
with increasing T to form the full Fermi surface at T* (right panel). From |4T| . 



conductivity measurements [40j . that yield similar asymmetries as in ARPES. 
Another generic feature of HTSC is a kink seen in the ARPES nodal dispersion 
as shown in figure[21Ia). This kink also effects the scattering rate of the coherent 
quasiparticles as measured by the momentum distribution curves (MDC) width, 
see figure IHb) and 

An interesting feature seen in ARPES is the shrinking of the Fermi surface 
when the pseudogap opens at T*. With decreasing temperature, more and 
more states around the antinodal region become gapped and the Fermi surface 
becomes continuously smaller. Instead of a full Fermi surface, the pseudogapped 
state exhibits Fermi arcs [321 [331 1331 STJ US] , that finally collapse to single nodal 
Fermi points at T = Tc (see figure fTO)) . For a detailed discussion on this and 
related ARPES observations, we refer the reader to the ARPES reviews in the 
literature p3l [24]. 

1.3.2 Scanning tunnelling microscopy (STM) 

In contrast to ARPES, STM is a momentum integrated probe. However, its 
ability to measure the local density of occupied as well as unoccupied states with 
a high energy resolution gives very valuable insights into HTSC. An example 
for a STM study of Bismuth-based HTSC is shown in figure [11] The data in 
the superconducting state reveals a density of states, which is characteristic of a 
d-wave gap, i.e. there is no full gap in contrast to s-wave superconductivity. In 
the pseudogap state (above Tc) the density of states is still suppressed around 
w = (zero voltage), however, the characteristic peaks disappear. Another 
interesting feature seen in figure [TT] is the striking asymmetry between positive 
and negative voltages, which becomes more pronounced for the underdoped 
sample. An explanation for this generic property of HTSC will be discussed in 
detail in the following sections. 

A key advantage of STM is the possibility to obtain spatial information. For 
example, STM experiments allow for the investigation of local electronic struc- 
ture around impurities [46l [47l [48] , and around vortex cores [49l [50l [51] in the 
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Figure 11: STM data for underdoped (UD) and overdoped (OD) Bi2212, and 
overdoped Bi2201; comparison between the pscudogap (dashed line, T > Tc) and 
the gap in the superconducting state (sohd hue, T < Tc). The underdoped data 
exhibit a significant asymmetry between positive and negative bias voltages. For 
an analysis of the temperature-dependent pscudogap see figure [T] From |36j . 
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superconducting state. Two other interesting features recently reported by STM 
are a checkerboard hke charge density wave [52l [53] and the existence of spatial 
variations in the superconducting gaps |54j . The origin of these observations is 
currently being debated intensely. 

1.4 Theories 

It is beyond the scope of this article to provide an overview of various theo- 
ries of high temperature superconductivity that have been put forward in the 
literature. Due to the enormous complexity of the experimentally observed fea- 
tures, it is not easy to agree on the key ingredients necessary for setting up a 
comprehensive theory. Further, the decision to trust new experimental results 
is often difficult, since the sample quality, experimental resolution, and the way 
the data is extracted arc often not completely clear. Not surprisingly perhaps, 
these circumstances have allowed for diverse theoretical approaches, motivated 
respectively by distinct aspects of the HTSC. In the following, we summarise a 
few theoretical approaches where the proximity of a superconducting phase to 
a Mott insulator and / or antiferromagnctism plays an important role. 

1.4.1 Electronic models 

To find an appropriate microscopic reference model is the first step in formulat- 
ing any theory. Such a model should be simple enough on the one hand to be 
treated adequately, but should also be complex enough to explain the relevant 
properties. In the case of the HTSC, it is widely accepted that strong correla- 
tions in the two-dimensional (2D) layers play an essential role. The copper- 
oxygen layers are appropriately described by a three-band Hubbard model, 
which includes the Cu da;2_j,2 -orbital and the two O p-orbitals [55l [56]. Its 
simplified version is an one-band Hubbard modeH, where each site corresponds 
to a copper orbital with repulsive on-site interaction between electrons |25j . The 
derivation of this model Hamiltonian can be found in the reviews of Lee, et al. 
and Dagotto [57]. 

1.4.2 Resonating valence bond picture 

Soon after the discovery of high Tc superconductivity, Anderson suggested the 
concept of a resonating valence bond (RVB) state [TS] as relevant for the HTSC. 
In this picture, the half-filled Hubbard model is a Mott insulator with one elec- 
tron per site. The charged states, doublons and holons, form bound charge- 
neutral excitations in the Mott insulating state and leads to the vanishing of 
electrical conductivity. Equivalently one can talk of virtual hopping causing a 
supcrcxchangc interaction J between the electrons at the copper sites. There- 
fore, the half-filled systems can be viewed as Heisenberg antiferromagnet with 
a coupling constant J. 

^Henceforth we refer to the one-band Hubbard model by the phrase "Hubbard model". 
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Anderson proposed that upon doping quantum fluctuations melt the antifer- 
romagnetic Neel lattice and yield a spin liquid ground state (denoted the RVB 
state) in which the magnetic singlet pairs of the insulator become the charged 
superconducting pairs. We will show in the following sections that the RVB 
picture provides a natural explanation for several key features of the HTSC 
such as the d-wave pairing symmetry, the shape of the superconducting dome, 
the existence of a pseudogap phase, the strong deviations from the BCS-ratio, 
and the singular k-dependence of the one-particle self-energy when approaching 
half-filling. 

1.4.3 Spin fluctuation models 

While the RVB idea approaches the problem from the strong coupling limit, i.e. 
large on-site electron repulsion U, spin fluctuation modcl^ start from the weak 
coupling (small U) limit. The technique extends the Hartrcc Fock (HF) random 
phase approximation and leads to a pairing state with d-wave symmetry. Within 
this picture, superconductivity is mediated by the exchange of antiferromagnetic 
spin fluctuations. 

Weak-coupling approaches such as spin fluctuation models essentially remain 
within the context of Landau theory of Fermi liquids for which the quasiparticle 
renormalisation is Z = m/m* , when the self-energy is not strongly k-dependent. 
Here, m* ~ and m is the bare band mass. The Fermi liquid relation Z ^ vp, 
is however difficult to reconcile with experimental results for the HTSC, as 
Z ^ and vp const for doping x — > 0, as we will discuss in more detail in 
section 16.1.41 

1.4.4 Inhomogeneity-induced pairing 

Within this class of theories, the proximity of high temperature superconduc- 
tivity to a Mott insulator plays an important role. It is postulated that the 
superconducting pairing is closely connected to a spontaneous tendency of the 
doped Mott insulator to phase-separate into hole rich and hole poor regions at 
low doping. The repulsive interaction could then lead to a form of local super- 
conductivity on certain mesoscale structures, "stripes" . Calculations show that 
the strength of the pairing tendency decreases as the size of the structures in- 
creases. The viewpoint of the theory is as follows: Below a critical temperature, 
the fluctuating mesoscale structures condense into a global phase-ordered su- 
perconducting state. Such a condensation is facilitated if the system were more 
homogeneous, however, more homogeneity leads to larger mesoscale structures, 
and thus weaker pairing. Therefore, the optimal Tc is obtained at an optimal 
inhomogcneity, where mesoscale structures are large enough to facilitate phase 
coherence, but also small enough to induce enough pairing. Within the phase- 
separation scenario spontaneous inhomogencities tend to increase even in clean 
systems when approaching half-flUing. In this framework, the pseudogap in the 

^For more details we refer to the review articles by Moriya and Ueda |58j . Yanase, et al. 
[591 . and Chubukov, et al. [60]. 
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underdoped regime can be understand as a phase, which is too granular to ob- 
tain phase coherence, but has strong local pairing surviving above T^. These 
ideas are reviewed in detail by Kivclson and collaborators [IHl EB El] • 

1.4.5 50(5) - theory 

Motivated by the vicinity of antiferromagnetism and superconductivity in the 
phase diagram of the HTSC, the 50(5)-theory [63j attempts to unify these 
collective states of matter by a symmetry principle. In the 5*0(5) picture, the 5 
stands for the five order parameters used to set up the theory; three degrees of 
freedom for antifcrromagnctic state {N^, Ny, N^) and two degrees of freedom 
for the superconducting state (real and imaginary parts of the superconducting 
order parameter). The theory aims to describe the phase diagram of HTSC 
with a single low-energy effective model. A so-called projected 50(5)-theory 
has been proposed to incorporate strong correlation effects. Several studies 
have also examined the microscopic basis for the 50(5) theory (see review by 
Demler, et al. [63]). 

1.4.6 Cluster methods 

Though numerical methods such as Lanczos (exact diagonalisation) and quan- 
tum Monte Carlo have been very popular [ST], they are limited by the (small) 
cluster size. All statements concerning the thermodynamic limit become impre- 
cise due to significant finite size effects. The "quantum cluster" method which 
aims to mitigate finite size effects in numerical methods, has been used by sev- 
eral groups to study strongly correlated electronic systems. These methods treat 
correlations within a single finite size cluster explicitly. Correlations at longer 
length scales are treated either perturbatively or within a mean field approxima- 
tion [64] • In recent years, this method has been used in several studies to extract 
the ground state properties of the Hubbard model. They reproduce several fea- 
tures of the Cuprate phase diagram and report d-wave pairing in the Hubbard 
model. However, even these sophisticated numerical methods are not accurate 
enough to determine the ground state of the Hubbard model unambiguously. 

1.4.7 Competing order 

In most of the theories outlined above, the pseudogap phase is characterised 
by the existence of preformed pairs. Hence, there are two relevant temperature 
scales in the underdoped regime. Pairs form at a (higher) temperature T*, and 
the onset of phase coherence at Tc leads to superconductivity. However, there 
are other theories that take the opposing point of view; viz., the pseudogap and 
superconductivity are two phases that compete with each other. In these sce- 
narios, the pseudogap is characterised by another order parameter, e.g. given 
by an orbital current state [65] or a d-density wave [66]. Thus, the pseudogap 
suppresses superconductivity in the underdoped regime, and can also partially 
survive in the superconducting state. These approaches predict that the pseu- 
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doping doping 



Figure 12: Two proposed theoretical phase diagrams for the Cuprates. a) RVB 
picture; b) Competing order scenario: the pseudogap (PG) ends in a quantum 
critical point (black dot); pseudogap and superconducting state (SC) can coexist 
(SC+PG). 

dogap line ends in a quantum critical point inside the superconducting dome. 
These two scenarios are contrasted in figure [121 

1.4.8 BCS-BEC crossover 

In this picture, the pseudogap is explained by a crossover from BCS to Bose 
Einstein condensation (BEG) [67l[68]- While in the BGS limit the fermionic 
electrons condensate to a superconducting pair state, the BEG limit describes 
the condensation of already existing pairs. In the crossover regime, one expects 
a behaviour very similar to that observed in the pseudogap of HTSG; formation 
of pairs with a corresponding excitation gap occurs at a temperature T* , and 
the pairs condense at a lower temperature Tc < T*. It is interesting to note that 
the physics behind this idea can be described by a generalisation of the BGS 
ground state wave function, |^'o) [68] • It is however unclear how to incorporate 
the antiferromagnetic Mott-Hubbard insulating state close to half-filling within 
a BGS-BEG crossover scenario. 

2 Resonating valence bond (RVB) theories 

The resonating valence bond (RVB) state describes a liquid of spin singlets, 
and was proposed originally as a variational ground state of the spin S ~ 1/2 
Heisenberg model (which describes the low energy physics of the Hubbard model 
at half filling). Anderson originally proposed that the magnetic singlets of the 
RVB liquid become mobile when the system is doped and form charged su- 
perconducting pairs. As we will discuss in this section, this idea has led to a 
consistent theoretical framework to describe superconductivity in the proximity 
of a Mott transition. In this section, we will discuss possible realisations of 
RVB superconductors along with the predictions of the theory. We also give 
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Figure 13: Hopping processes with a virtual doubly occupied site corresponding 
to the SfSj and S" term of the Heiscnberg Hamiltonian. respectively; in the 
case of parallel spins virtual hopping is not possible. 

a outlook on the implementations of the RVB picture by Gutzwiller projected 
wave functions, slave boson mean field theory and the bosonic RVB approach. 

2.1 The RVB state - basic ideas 

Within the resonating valence bond (RVB) picture, strong electron correlations 
are essential for superconductivity in the Cuprates. The Hubbard model is 
viewed as an appropriate microscopic basis and the corresponding many-body 
Hamiltonian is given by, 

H = - hv) {4aC3<y + c]aC^i<y) + ^ X! "*T"U : (1) 

where c\ creates and q annihilates an electron on site i. The hopping integrals, 
connect sites i and j. We shall restrict our attention to nearest neighbour 
hopping t for the moment and will also discuss the influence of additional hop- 
ping terms subsequently. The operator riio- = c\^Cic, denotes the local density 
of spin a =4, ] on site i. We consider an on-site repulsion U ^ t, i.e. we work 
in the strong coupling limit, which is a reasonable assumption for the HTSC. 

2.1.1 RVB states in half-filled Mott-Hubbard insulators 

Let us first consider the half-filled case. Since U is much larger than t the mean 
site occupancy is close to charge neutrality, namely one. It costs energy U for 
an electron to hop to a neighbouring site. This potential energy is much higher 
than the energy the electron can gain by the kinetic process. Thus, the motion of 
electrons is frozen and the half-filled lattice becomes a Mott-Hubbard insulator. 
However, there are virtual hopping processes, where an electron hops to its 
neighbouring site, builds a virtual doubly occupied site, and hops back to the 
empty site. Such virtual hoppings lower the energy by an amount of the order 
J = At^ /U . Pauli exclusion principle allows double occupancy only for electrons 
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Figure 14: Left: Antifcrromagnetic Neel lattice with some holes. The motion 
of a hole (consider bold circles) frustrates the antifcrromagnetic order of the 
lattice. Right: Snapshot of the RVB state. A configuration of singlet pairs 
with some holes is shown. The RVB liquid is a linear superposition of such 
configurations. 

with opposite spin (see figure [13]). Thus, virtual hopping favours anti-parallel 
spins of neighbouring electrons and we obtain an effective antifcrromagnetic 
Heiscnbcrg Hamiltonian, 

= J ^ Si • Sj , J > , (2) 

with an antifcrromagnetic exchange constant J = At'^/U, the spin-operator Si 
on site i, and (ij) denoting a sum over nearest neighbour sites. At the level 
of mean field theory, i.e. treating the spins semiclassically, the 2D Heiscnbcrg 
model on a square lattice has an antifcrromagnetic Nccl ground state with long 
range order and broken symmetry (left panel of figure I14p . This molecular- 
field prediction is experimentally (by neutron scattering studies [69]) as well as 
theoretically (by a quantum nonlinear a model |70| ) well established. 

Anderson [fs] suggested that a resonating valence bond (RVB) liquicill is 
very close in energy to the Neel state for undoped Cuprates. Instead of a Neel 
state with broken symmetry, a fluid of singlet pairs is proposed as the ground 
state; i.e., the ground state is described by a phase coherent superposition of 
all possible spin singlet configurations (see right panel of figure [T4]) . For spin 
S = 1/2, quantum fluctuations favour such singlets than classical spins with 
Neel order. To see this, consider a one-dimcnsional (ID) chain (sec figure [TS]) . 
In this case, a Neel state with Sz ~ ±1/2 gives an energy of — J/4 per site. 
On the other hand, the ground state of two antifcrromagnetic coupled spins 
5* = 1/2 is a spin singlet with — S {S + 1) J = —3 / 4 J . It follows that a chain of 
singlets (see figure [T5]) has an energy of —3/8 J per site, much better than the 
Nccl-ordcred state. This simple variational argument shows that a singlet state 
is superior in ID. Similar considerations for the 2D Heiscnbcrg model give the 

^Long before the discovery of HTSC Anderson and Fazckas 1131 1141 proposed the RVB 
hquid as a possible ground state for the Heisenberg model on a 2D triangular lattice. 
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Figure 15: Neel state (left) and singlet state (right) for ID antiferromagnetic 
spin S — 1/2 chain. 



energies —1/2 J per site for the Neel lattice, the singlet state remains at —3/8 J 
per site. Following this reasoning we find that singlets become much worse than 
the Neel state in higher dimensions. 

Liang, Doucot, and Anderson [71] showed that the singlet 'valence bonds' 
regain some of the lost antiferromagnetic exchange energy by resonating among 
many different singlet configurations and become therefore competitive with the 
Neel state in 2D. The resonating singlets are very similar to benzene rings with 
its fluctuating C-C links between a single and a double bond; an analogy that 
motivated the term 'RVB'. 

2.1.2 RVB spin liquid at finite doping 

Though an antiferromagnetically long range ordered state is realised in the 
undoped insulator, the order melts with only a few percent of doped holes. To 
understand this, consider the example shown in figure [T3] (left). The figure shows 
that moving holes cause frustration in the antiferromagnetic but not the RVB 
state, figure fT4l (right). A single hole moving in the background of a Neel state 
was studied extensively by several author^, and analytical calculations showed 
that the coherent hole motion is strongly renormalised by the interactions with 
the spin excitations [751 [73]. When more holes are injected into the system, 
the interaction of the holes with the spin background completely destroys the 
antiferromagnetic Neel state and an RVB liquid (or spin liquid) state becomes 
superior in energy. Then the singlet pairs of the RVB liquid are charged and 
may condense to a superconducting ground state. 

2.2 Realisations and instabilities of the RVB state 

Whether there exist two dimensional models with an RVB ground state is still 
an open question. We may however regard the RVB state as an unstable fix- 
point [74] prone to various instabilities. The situation is then analogous to that 
of the Fermi liquid, which becomes generically unstable in the low-temperature 
limit either towards superconductivity or various magnetic orderings. For in- 
stance, Lee and Feng studied numerically how a paramagnetic RVB state can 

®The single hole problem together with the corresponding literature is discussed in |19) in 
more detail. 
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be modified to become a long range (antiferromagnetically) ordered state by 
introducing an additional variational parameter [75j . In this view of antiferro- 
magnetism, the "pseudo Fermi surface" of the insulating RVB state undergoes a 
nesting instability to yield long range antiferromagnctic order |761 177] . In figure 
[TBI we present an illustration of the concept of the RVB state as an unstable 
fixed point. In the following, we discuss this point further. 

Besides the square lattice with nearest neighbour hopping, the RVB spin 
liquid was proposed as a ground state on a square lattice with further neigh- 
bour hopping as well as in a triangular lattice. Experiments [7S] indicate that 
such a spin liquid state may be realized in the organic compound k-(BEDT- 
TTF)2Cu2(CN)3 which is an insulator in the proximity of a Mott transition. 
Trial spin liquid wave functions using Gutzwiller projected RVB states have 
been proposed in this context by Motrunich [79]. A U{1) gauge theory of the 
Hubbard model has also been invoked to study this system [80]. Although the 
simple Neel ordered state is destroyed due to frustration in these cases, the 
RVB spin liquid (at n = 1) does not become the (T = 0) ground state, which is 
either a valence bond crystal state [STJ [511 IMl IHll [5S] or a coplanar 120° anti- 
ferromagnetic ordered state [86], respectively. In addition, instabilities against 
inhomogeneous states like stripes [IHl IMl IHl] are conceivable, and are not ex- 
phcitly included in figure [TBI A recent ARPES study on La2_2,Ba3,Cu04 [57], 
which exhibits static charge order and suppressed superconductivity around 
doping X = 1/8, supports the idea that the superconducting RVB state can be 
continuously connected and unstable against a charge ordered state. 

Nevertheless an RVB state can be realised if a finite number of hole is in- 
duced into the system, viz. when the bosonic spin state realised at half-filling 
turns into a free fermionic state by the introduction of charge carriers. The 
hopping processes then destroy above instabilities towards magnetic or valence 
bond crystal ordering and a superconducting RVB state can be stabilised. A 
schematic picture of this scenario is presented in figure 1161 

In the case of HTSC, holes are created by changing the doping concentra- 
tion. A similar mechanism was proposed for superconductivity in the triangular 
lattice based Cobaltates [55] [53] . Within RMFT calculations such a triangular 
model would result a d -I- jc?-wave pairing state [HO]. On the other hand, an RVB 
superconducting state at half-filling just below the Mott transition [91] was re- 
cently suggested for organic superconductors [92l [93] [94| ■ Here, the necessary 
holes could result from a finite number of conducting doubly occupied sites as 
illustrated in figure [16l 

To summarise, an RVB superconductor could emerge by two different mech- 
anisms starting from a Mott insulating system {n — 1 and U > Uc); cither 
upon doping (n ^ 1), or from self doping a half-filled system close to the the 
Mott-Hubbard transition [U ~ Uc)- In this review, we focus our attention on 
the former possibility, i.e. the occurrence of an RVB superconductor in a doped 
Mott-Hubbard insulator. 
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Figure 16: Schematic picture of instabilities and realisations of the RVB spin 
liquid state, viz. of the RVB state as an unstable fixpoint. The top panel 
shows an RVB spin liquid at half-filling in the Mott-Hubbard insulating limit 
(C/ > Uc)- The middle panel illustrates instabilities of the RVB liquid state 
in a square lattice, a frustrated square lattice, and a triangular lattice in the 
half-filled limit. The lower panel shows realisations of the RVB liquid, which 
are realised at finite doping or close to the Mott-Hubbard transition ([/ ~ U^. 
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doping 

Figure 17: RVB phase diagram with singlet pairing temperature T* and phase 
coherence temperature Tcoh (superconducting state SC, pseudogap PG). 

2.3 Predictions of the RVB hypothesis for HTSC 

In this subsection wc discuss some predictions from RVB theory, which agree 
weh with experimental observations. As wc will show in the following sections, 
the arguments we present here are substantiated by more detailed microscopic 
calculations. 

Within the RVB picture, a possible explanation for the temperature-doping 
phase diagram is obtained by considering two temperature scales (figure [T7|) . 
The singlets of the RVB liquid form at temperature T* , a temperature scale 
which decreases away from half-filling [95| owing to the presence of doped and 
mobile holes. Holes, on the other hand, allow for particle number fluctuations, 
which are fully suppressed at half-filling, and thus enhance the stability of the 
superconducting state against thermal fluctuations. This results in a second 
temperature, Tcoh, which increases with doping and below which the supercon- 
ducting carriers become phase coherent. The superconducting transition tem- 
perature Tc is therefore determined by the minimum of T* and Tcoh as shown 
in figure [TTlinS]. 

It is evident from the above picture that a pseudogap forms for Tcoh < T < 
T*, i.e. for underdoped samples. In this state, although phase coherence is lost, 
the RVB singlet pairs still exist. Therefore, we have to break a pair to remove 
an electron from the copper-oxygen layers within the pseudogap regime. The 
resulting excitation gap manifests itself, e.g. in the c-axis conductivity or in 
ARPES measurements. 
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These schematic explanations are confirmed to a certain extent by analyti- 
cal as well as numerical calculations (at zero temperature). RMFT and VMC 
methods show an increase of the superconducting gap, but a vanishing super- 
conducting order parameter, when approaching half-filling. This behaviour is 
in complete agreement with the T observations in experiments. It also 
explains the strong deviation from the BCS-ratio in the undcrdoped regime of 
the HTSC, if the superconducting order parameter is related to Tc- On the 
other hand, the doping dependence of the onset temperature of the pseudogap 
T* can be related to the magnitude of the gap at T = (in agreement with 
experiments, see figure [T]). 

Perhaps the most remarkable prediction of the RVB theory was the c?-wave 
nature of the superconducting state. A d-wave superconducting state was pre- 
dicted by RVB based studies as early as in 1988 [H [96l [17l [97] , long before 
the pairing symmetry was experimentally established. These early calculations 
also correctly described the vanishing of superconductivity above about 30% 
doping. 

Implementing the RVB idea by projected wave functions, one finds a nat- 
ural explanation of the suppression of the Drude weight and of the superfluid 
density in the underdoped regime as well as the particle-hole asymmetry in 
the density of single particle states. Further successes of the RVB theory are 
calculations that predict a weakly doping dependent nodal Fermi velocity, but 
a quasiparticle weight which is strongly doping dependent- the quasiparticle 
weight decreases with doping x in agreement with ARPES experiments. These 
effects can be understood by a decrease in the density of freely moving carriers 
at low doping, which results in a dispersion mainly determined by virtual hop- 
ping processes (proportional to the superexchange J). In the half-filled limit, 
this behaviour results in a divergence of the k-dependence of the electron's self- 
energy, lim^^o f^5](w, k = kp')/9<x' ~ 1/x oo , which transcends the nature 
of orthodox Fermi liquids. These will be discussed in more detail in sections [6] 
and El 

In addition to the above key features of HTSC, RVB theory has also been 
successfully applied to several other phenomena such as charge density patterns 
[Ml [Ml llOOl I101| , the interplay between superconductivity and antiferromag- 
netism [lOl [Ml [IM [IM [M], impurity problems [TOl [Ml [TTO], and 
vortex cores |11H I112| . 

In conclusion, analytical and numerical results provide significant support 
to the RVB concept. However, most RVB studies are restricted to zero temper- 
aturc@, making the finite temperature picture detailed above, somewhat spec- 
ulative. Extending the calculations to finite temperature is an important and 
open problem in the theory of RVB superconductivity. A related issue is the 
destruction of superconductivity in the underdoped samples where we expect 
phase fluctuations to play an increasingly important role at low temperatures 
[1141 1115] since particle number fluctuations are frozen in the proximity of the 

possible ansatz for finite temperatures was recently proposed by Anderson 11131 . He 
suggests a spin-charge locking mechanism within the Gutzwiller-RVB theory to describe the 
pseudogap phase in the underdoped Cupratcs as a vortex liquid state. 
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Mott insulator. It is presently an unsettled question to which extent this picture 
is equivalent to alternative formulations, such as an increase of inhomogeneities 
(like in the 'inhomogcneity- induced pairing' picture [551 [Ml El] ) or a destruc- 
tion of the superfluid density due to nodal quasiparticle excitations (see section 
16. 4p . that were also proposed to describe the transition from the superconduct- 
ing state to the pseudogap state in the undcrdopcd regime. Further work is 
necessary to clarify this point. 



2.4 Transformation from the Hubbard to the t-J model 

The RVB scenario is based on the existence of a strong antifcrromagnetic su- 
perexchangc, J. The superexchangc process via virtual hopping processes re- 
sults in an effective Hciscnbcrg Hamiltonian as discussed earlier (see figure I13p . 
We now present a more formal and systematic derivation of a low energy theory 
starting from the Hubbard Hamiltonian in the strong coupling limit {U ^ t). 
The basic idea is to make the theory "block diagonal" ; i.e., subdivide the Hamil- 
tonian matrix elements into processes that preserve the local number (diagonal 
processes) and those that do not (off-diagonal) by a unitary transform. Since 
we will be interested in the strong coupling limit, off-diagonal processes will 
be removed as such (high energy) configurations are not allowed in the Hilbert 
space of the effective (low energy) theory. 

The unitary transformation, e"*"^ to lowest order in t/U [1161 1117] can be 
obtained as follows. First we assume that S is of the order 0{t/U) and expand 
the transformed Hamiltonian, 

^ e^'^He-''^ ^e:^{f + U)e-'^ (3a) 

= T + U + i[S, T + U] + - [S,[S,T+U]] + ... (3b) 

= U + T + t[S,U]+i[S,T] + - [S,[S,U]]+ ^ . (3c) 

Here, we split the Hubbard Hamiltonian H into the kinetic energy part T, 
the first term of ([T]), and the potential energy part [/, the second term of ([T]) 
(includes the parameter U). In ([3c[) we have ordered the terms in powers of t/U. 
For a block diagonal Hamiltonian to order 0{t/U), the term, T + i [S, U], 

in ([5c[) may not contain any (real) hopping processes changing the total number 
of doubly occupied sites. An appropriate choice for S is given by. 



S 



E ^-Jr{<^Ld,,<r + ''Ld^,. - h.c.) , (4) 
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f + z[S,U]^-Y. [ala,, + did,, + h.c. j , (5) 
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does not involve hopping process changing the number of double occupancies. 
Here, we used the operators a| ^ = (1 — ni__cr)cj ^ and dj ^ = ni_-ac\^. Equation 
([5]) is block diagonal and verifies the choice of 5 in 

The full form of is now obtained by evaluating all C(t^/?7)-terms 

in (Pc|) with S from By restricting ourselves to the subspace of no double 
occupancies (the low energy subspace or the lower Hubbard band), wc find the 
t-J Hamiltonian, 

Ht-j = Pg H^'f^^ PG = PG{T + Hj + Hi)PG , (6) 

where, 

PG-^(l-n,Tn,x) , (7) 

i 

is the Gutzwiller projection operator that projects out all doubly occupied sites. 
The terms of the Hamiltonian are given by, 

^ = - {4,aCj,a + S>Ci,<T) , (8) 

Hj = ^ J(,j) fs, - in,nj j , (9) 

J(i+Ti,iA+T2) t t 



^3 = - E 



E 



J(i+Ti,i.i+T2) t t fl(\\ 



where J{ij) = ^'tfij)/U and J{i,j,i) = '^t{i,j)t(j.i)/U . arc pairs of neighbour 

sites and « + T(i 2) denotes a neighbour site of i. Equation ([6]), together with ([8])- 
(fTU|) . gives the full form of t-J Hamiltonian. However, the so-called correlated 
hopping or three-site term is often ignored since its expectation value is 
proportional both to t'^ /U and the doping level x. Further, the density-density 
contribution riirij is sometimes neglected within the superexchange term Hj, 
as it is a constant at half-filling. Note that ([8]) is equivalent to (O due to the 
projection operators Pg occurring in the definition ([6|) of the t-J Hamiltonian. 

The unitary transformation illustrates the relationship between superex- 
change and the physics of the (strong coupling) Hubbard model. We see that 
as a result of the unitary transform, the low energy model is given by the t-J 
Hamiltonian ([6]) which does not allow for double occupancies. At half filling, 
each site is singly occupied and the hopping of electrons is frozen since real hop- 
ping now leads to states in the upper Hubbard band. As a result, the kinetic 
energy term in the Hamiltonian vanishes, and the t-J Hamiltonian reduces to 
an antiferromagnetic Heisenberg model 

The original Hamiltonian relevant for the cuprates contains three bands 
per unit cell, one copper band and two oxygen-derived bands. One band only 
crosses the Fermi surface with a single effective degree of freedom per unit cell. 
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the Zhang-Rice singlet [21], corresponding to an empty site in J terminology. 
Using this venue, the hopping matrix elements and the superexchange param- 
eters relevant for the t — J model could be derived directly. The Hubbard-f/ 
entering the relations derived above then takes the role of an effective modelling 
parameter. 

2.5 Implementations of the RVB concept 

The t-J Hamiltonian ([6]) is more suitable than the Hubbard model for studying 
RVB superconductivity, because it includes the superexchange term explicitly, 
and it is this term which is responsible for the formation of singlets. However, 
for exact numerical methods, the t-J Hamiltonian provides only a minor sim- 
plification over the Hubbard Hamiltonian, and one must turn to approximate 
schemes for any calculations on sufficiently large clusters. In the following, we 
start with the t-J Hamiltonian as an appropriate microscopic model for HTSC, 
and briefly discuss three schemes that allow for systematic calculations of the 
RVB state. 

2.5.1 Gutzwiller projected wave functions 

Anderson proposed projected BCS wave functions as possible RVB trial 
states for the t-J model. These states provide a suggestive way to describe an 
RVB liquid in an elegant and compact fornFl. 

|*rvb) =P7vPg|BCS) , (11) 

with the BCS wave function 

|BCS)=n("k + i^k4^cl^i)|0) , (12) 

k 

which constitutes a singlet pairing state. Here, the operator Pq (Gutzwiller pro- 
jection operator) projects out double occupancies and the P/v fixes the particle 
number to A^; Uk and Vk are the variational parameters with the constraint, 
^k + ^k = 1- The form of I^'rvb) provides a unified description of the Mott 
insulating phase and the doped conductor. It immediately suggests the pres- 
ence of singlet correlations in the undoped correlations and relates them to a 
superconducting state away from half filling. 

Projected wave functions were originally proposed by Gutzwiller in 1963 to 
study the effect of correlations presumed to induce ferromagnetism in transition 
metal compounds [118j . In subsequent years, these wave functions were applied 
to study the Mott-Hubbard metal insulator transition [119j and for a description 
of liquid ■^He as an almost localised Fermi liquid |116[ 11201 1121] . etc. However, 
these early studies considered only a projected Fermi sea, 

Pg|*fs) = n ^T^kilO) . (13) 

k<kF 

^"For a real space representation of equation we refer to section TS.!.!! 
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in the Hubbard model, whereas Anderson [TS] suggested a projected BCS paired 
wave function for the t-J model. 

To calculate the variational energy of a projected state l^*) = Pg\^o), ex- 
pectation values of the form 



{■^o\PGOPG\^n 
{^o\PgPg\'^q) 



(14) 



must be considered, where O is the appropriate operator. Here, j^'g) can be 
any wave function with no restriction in the number of double occupancies, viz. 
it lives in the so-called 'pre-projected' space. The choice of In our case we 
concentrate on |5'o) = |BCS). In section [Z6l we will review a few other types 
of trial wave functions used to study correlated electron systems. The exact 
evaluation of (fT4|) is quite sophisticated and requires variational Monte Carlo 
(VMC) techniques that will be discussed in section [5l However, approximate 
analytical calculations can be done by a rcnormalisation scheme based on the 
Gutzwillcr approximation (GA). The GA will will be outlined in the sections [3] 
andlH Within this approximation, the effects of projection on the state |\E'o) are 
approximated by a classical statistical weight factor multiplying the expectation 
value with the unprojcctcd wave function |120j . i.e. 

' 9o i.T, X • (15) 



{^o\PgPg\^o) (^oI^o) 

The so-called Gutzwillcr rcnormalisation factor go only depends on the local 
densities and is derived by Hilbcrt space counting arguments [20l I120[ 1122] or 
by considering the limit of infinite dimensions (d = oo) [123[ 11241 11251 [T26] . The 
GA shows good agreement with VMC results (see [20] ) and is discussed detailed 
in section [3l 

Gutzwillcr projected wave functions thus have the advantage that they can 
be studied both analytically (using the GA and extensions thereof) and nu- 
merically (using VMC techniques and exact diagonalisation) . Since these wave 
functions provide a simple way to study correlations such as pairing correla- 
tions, magnetic correlations etc. in the presence of a large Hubbard repulsive 
interaction, they have been used extensively in the literature. As we will show in 
the following sections, the Gutzwillcr- RVB theory of superconductivity explains 
several key features of the HTSC . More generally, we believe this approach 
is sufficiently broad that it could be used to study a wide range of physical 
phenomena in the proximity of a Mott transition. 



2.5.2 Slave boson mean field theory (SBMFT) and RVB gauge the- 
ories 

Another representation of the t-J Hamiltonian, equation is obtained by re- 
moving the projection operators Pgi a-nd replacing the creation and annihilation 
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operators by 

4<T ^ c\.A.a (1 - n.,,-<y) ' and (16a) 

Ci,cr Ci^cr = Ci^a (1 - ?^^,-(T) , (16b) 

with (T =1, 4 and — cr denoting the opposite spin of <t. In this form the restriction 
to no double occupation is fulfilled by the projected operators c\ ^ and Ci^n- 
Thus, only empty and single occupied sites are possible, which can be expressed 
by the local inequality 

cr 

However, the new operators do not satisfy the fermion commutation relations, 
which makes an analytical treatment difficult. The slave-boson method [127[ 
112811129] handles this problem by decomposing c\ ^ into a fermion operator fl ^ 
and a boson operator hi via 

4 = flah ■ (18) 

The physical meaning of fj^ (fi.a) is to create (annihilate) a single occupied 
site with spin ct, those of bi (fo|) to annihilate (create) an empty site. Since 
every site can either by single occupied by an "f-electron, single occupied by a 
J,-electron, or empty the new operators must fulfill the condition 

fy^1+Anl+blb^ = l ■ (19) 

When writing the Hamiltonian in terms of the slave fermion and boson operators 
the constraint p9p is implemented by a Lagrangian multiplier A^. In the slave- 
boson representation, the t-J model is thus written as, 

Ht-j = - J2 {flMh- + flMk^r) (20) 

where the Heisenberg exchange term, 

is a function of fermion operators only, since superexchange does not lead to 

charge fluctuations. [95]. Furthermore, a chemical potential term, —/io X^i o- flafi.o-^ 
is included within the grand canonical ensemble. 

The advantage of this representation is that the operators {fia, bi) obey 
standard algebra and can thus be treated using field theoretical methods. The 
partition function Z of (|20p can be written as a functional integral over coherent 



30 



Bose and Fermi fields, allowing to calculate observables in the original Hilbert 
space. The Fermi fields can be integrated out using standard Grassmann vari- 
ables. Then carrying out a saddle-point approximation for the Bosc fields repro- 
duces the mean field level. The incorporation of Gaussian fluctuations around 
the saddle point approximation provides a possibility for systematic extensions 
of the SBMFT. One way to implement the constraint of single occupancy is to 
formulate the problem as a gauge theory. 

The development of RVB correlations and a superconducting phase in a 
lattice model as a gauge theory was first studied by Baskaran and Anderson [130j 
These authors noted that the Heisenberg Hamiltonian has a local U{1) gauge 
symmetry, which arises precisely because of the constraint of single occupancy. 
One may then develop an effective action which obeys this local symmetry and 
use it to calculate various averages. Since the free energy exhibits the underlying 
gauge symmetry, it is possible to go beyond mean field theory when calculating 
averages of physical quantities. Doping turns the local gauge symmetry into 
a (weaker) global U{1) symmetry which can be broken spontaneously, leading 
to superconductivity. Subsequently, Wen and Lee introduced an SU{2) gauge 
theory which leads to RVB correlations and superconductivity in a doped Mott 
insulator |131| . These approaches are reviewed in a recent work by Lee, et 
al. [H]. It should be noted that the Gutzwiller approximation, the SBMFT 
(which is the mean field solution about which gauge theories are constructed) 
are similar in the sense that both model the doped Mott insulator. In particular, 
real kinetic energy is frozen as one approaches half filling, and enhanced RVB 
correlations. In general, the results from SBMFT are quite similar to those 
from RMFT, e.g. the early prediction of d-wave superconductivity in the t-J 
model rests on very similar gap equations in both schemes. The SBMFT result 
showing d-wave pairing by Kotliar and Liu [TS] and by Suzumura, et al. [SS] 
nearly simultaneously appeared with the respective RMFT study by Zhang, et 
al. [20j. These studies followed an earlier work of Baskaran et al. [95], who 
initially developed a slave boson theory for the t-J model. For a more detailed 
review on SBMFT we refer to The SBMFT and Gutzwiller approaches 

differ in the way the local constraint is treated and consequently, there are 
quantitative discrepancies between these approaches. Some of these will be 
highlighted in subsequent sections of this review. 

2.5.3 The bosonic RVB theory 

As the name indicates, this approach is based on a bosonic description of the 
t-J model. The advantage of this method is that it accounts well for the an- 
tiferromagnetic correlations of the Heisenberg model at half filling as well as 
of the hole doped t-J model. At half filling, the ground state of the bosonic 
RVB theory (b-RVB) is related to the RVB wave function of Liang et al. [7T] 
which is the best variational wave function available for the Heisenberg model. 
The basic premise of the b-RVB theory is that hole doping an insulator with 
AF correlations (not necessarily long ranged) lead to a singular effect called the 
"phase string" effect }132j . A hole moving slowly in a closed path acquires a 
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nontrivial Berry's phase. As this effect is singular at the length scales of a lattice 
constant, its topological effect can be lost in conventional mean field theories. 
So, the theory proposes to take this effect into account explicitly before invoking 
mean field like approximations. The electron operator is expressed in terms of 
bosonic spinon and holon operators, and a topological vortex operator as. 

The phase operator Oio- is the most important ingredient of the theory and re- 
flects the topological effect of adding a hole to an AF background. The effective 
theory is described by holons and spinons coupled to each other by link fields. 

Away from half filling, the ground state of the b-RVB theory is described 
by a holon condensate and an RVB paired state of spinons. The supercon- 
ducting order parameter is characterised by phase vortices that describe spinon 
excitations and the superconducting transition occurs as a binding/unbinding 
transition of such vortices [133| . The theory leads naturally to a vortex state 
above Tc of such spinon vortices |134| . Bare spinon and holon states are confined 
in the superconducting state and nodal (fermionic) quasiparticles are obtained 
as composite objects Ll3_5,. 

The b-RVB theory realises transparently, the original idea of Anderson of 
holes moving in a prepaired RVB state. As mentioned above, the theory leads 
to definite and verifiable consequences such as a vortex state of spinons above Tc 
and spinon excitations trapped in vortex cores. However, the exact relationship 
between the b-RVB ground state and the simple Gutzwiller projected BCS wave 
function has not been clarified yet [136] . 

2.6 Variational approaches to correlated electron systems 

In this subsection, we briefly discuss how projected states, 

l^-) - Pg|*o) , (21) 

can be extended to study a wide variety of strongly correlated systems, varia- 
tional basis. Apart from the HTSC, these wave functions have been used in the 
description of Mott insulators |137| . superconductivity in organic compounds 
[MlIISHj, Luttingcr liquid behaviour in the t-J model [l39lfT40] . 

2.6.1 Order parameters 

A simple extension of the trial state pT|) is to allow for additional order pa- 
rameters in the mean field wave function \'^q). In section B.5.1[ we restricted 
ourselves to a superconducting BCS wave function |^'o) = |BCS). However, 
antiferromagnetic [75l [Ml [1021 [lOll [105] , ^r-flux [MTI [Ml [112! , or charge or- 
dered [98 1 (99 1 [Tool I lOlj niean field wave functions can also be used for |^'o). In 
addition, a combination of different kind of orders is possible. As an example, 
consider the trial wave function, 

|^o)=n("k + «kfo[T^Lki)|0) , (22) 

k 
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with 

bka = akCkcr + Cr/3kCk+Qcr • (23) 

Equation ()22|) includes finite superconducting as well as antiferromagnetic order 
|103j . Here, 6k is the Hartree-Fock spin- wave destruction operator with Q = 
(tt, tt) as required for a commensurate antiferromagnet. The parameters, and 
/3k are related to the antiferromagnetic order parameter Aaf by usual mean field 
relations; similarly, the superconducting order parameter determines the values 
of t'k and Uk- In the sections 2] and [SJ we will discuss applications of above wave 
function for the HTSC. 

We note that j'l'o) is applicable to all lattice geometries. It has been used, 
for instance, to study superconductivity in triangular lattice based Cobaltates 
[Ml [Ml [Ml [90] and organic compounds [92l [9l [Ml [IM]- Recent calculations 
show that projected states also provide a competitive energy on more exotic 
models such as a spin- 1/2 Heisenberg model on a Kagome lattice |143| . 

2.6.2 Jastrow correlators 

The incorporation of Jastrow correlator J7 jl44| provides an additional powerful 
way to extend the class of (projected) trial wave functions. In ([2T|) . the orig- 
inal Gutzwiller projector Pq can be viewed as the simplest form of a Jastrow 
correlator, 

Pc^Jg^ 5^' ""^"-'^ = n (1 - (1 - 9)n^,^^^,l) ■ (24) 

i 

So far we only considered Pq in the fully projected limit, which corresponds to 
g —f in Jg. However, when using (j24p in the Hubbard model, g becomes a 
variational parameter that determines the number of doubly occupied sites. 

The variational freedom of the trial wave function can be increased by in- 
cluding further Jastrow correlators, 

I*) = J. Jhd Jd Pg |*o) = Js Jhd Jd Jg |*o) ■ (25) 
Popular choices of Jastrow correlators arc the density-density correlator jTij, 

Jd = exp - ^ (1 - ni){l - Uj) , (26) 
the holon-doublon correlator Jhdi 

Jhd = exp - ^ Wij{h.idj + dihj) , (27) 
with hi — {\ — ni^){l — Tiii) and di ~ ni^riii, and the spin-spin correlator jTg, 

J, = exp[-^Uy-5f5H . (28) 
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The corresponding variational parameter are given by Vij , Wij , and Uij , respec- 
tively. 

Since the generalised trial wave function (|25p includes a very high number 
of variational parameters, one invariably chooses a small set depending on the 
problem at hand. In the case of the t-J model the situation is slightly simplified, 
because double occupancies are forbidden and thus 5^0 and Wij = 0. 

We now discuss the properties of the density-density correlator in ((26|) and 
assume Uij = Wij = for a moment. A positive Vij implies density-density repul- 
sion, a negative Vij means attraction and may lead to phase separation. Several 
studies indicate the importance of long range density-density Jastrow correlators 
for improving the variational energy. Hellberg and Mele [139] showed that the 
one-dimensional t-J model can be accurately described when ~ log \i — j\, 
i.e. when the Jastrow correlator is scale invariant. The incorporation of long- 
ranged density-density correlations induces Luttinger liquid like behaviour in 
the t-J model jl391 1140] . In the one-dimensional Hubbard model an appropri- 
ate choice of the density-density correlator in momentum space allows to distin- 
guish between metallic and insulating behaviour |137j . In the two-dimensional 
t-J model, Jd is often used to improves the variational energy of a projected 
superconducting state [1451 1146| as we will discuss in section 15.21 

The holon-doublon Jastrow correlator jT^fi is important for studying the re- 
pulsive Hubbard model on a variational basis. A negative Wi,j < implies 
attraction of empty and doubly occupied sites which ultimately may lead to 
a Mott-Hubbard insulating state (the Mott transition) [94l I138j . In two di- 
mensions, a negative nearest neighbour Wij ~ —Si^ij), substantially decreases 
the variational energy [Ml I138j , since these states occur as intermediate states 
during the superexchange process (compare figure fTS]) . Combining these effect 
with a superconducting wave function = |BCS) then explains key aspects 
of superconductivity in organic compounds near the Mott-Hubbard transition 
[94l 1138] . The Wij seems to be less important for one dimension, and is likely to 
be a consequence of the very good spin-spin correlation energy of the Gutzwiller 
wave function in one dimension [1161 1147] . 

The spin-spin Jastrow correlator J7s is not as often used as the density- 
density and the holon-doublon Jastrow correlators {Jd and Jhd)- However, 
recent studies show that the inclusion of is important when considering 
charge fluctuations within the two-dimensional t-J model |148j . An appropriate 
spin-spin Jastrow correlator J7s can also create antiferromagnetic order in a 
non-magnetic wave function, an example for the ability of Jastrow correlators 
to induce a new long-range order not manifest in the unprojected wave function. 

3 Gutzwiller approximation 

The Gutzwiller approximation (GA) is a straightforward method to handle 
Gutzwiller projected wave functions, that incorporate strong electron corre- 
lations by prohibiting doubly occupied sites. Within the GA, the effects of pro- 
jection are absorbed by statistical weight factors (Gutzwiller renormalisation 



34 



factors), which then allow for an analytical treatment of strongly correlated 
Gutzwiller wave functions. 

In this section, we present the derivation of the Gutzwiller factors by Hilbert 
space counting argument as well as considering the limit of infinite dimensions. 
Further, we discuss the importance of fugacity factors in the GA when com- 
paring analytical results with variational Monte Carlo (VMC) calculations in 
the canonical and grand canonical scheme, respectively. As we will show in the 
last part of this section, the GA can also be extended to the case of partially 
projected wave functions, where the projection operator does not act on a single 
"reservoir" site in the system. 

3.1 Basic principles of the Gutzwiller approximatiorJ"1 

The Gutzwiller approximation (or Gutzwiller renormalisation scheme) consti- 
tutes the basis of the RMFT and is a successful method to treat Hilbert space re- 
strictions due to strong electron correlations. It was applied originally |118illl9] 
to calculate the variational energy of the projected Fermi sea, Pg\FS), in Hub- 
bard like models. In these and other early papers, the projection operator, 
Pg = Yiii^ ~ c^^^iT^ii)' generalised to partial projection with the param- 
eter a determined by optimising the energy. Partially projected states were 
used successfully in modelling normal liquid ^He jl20[ I149| and heavy fermion 
systems [TCTlfTn] . 

Here, we focus on the t-J model {i.e., the large U limit of the Hubbard 
model). Consequently, we shall mainly discuss the fully projected case, i.e. 
a = 1. We will derive the corresponding renormalisation factors (Gutzwiller 
renormalisation factors) in this limit, and will not not discuss the generalisation 
to finite double occupancy. The latter case is obtained easily following the same 
reasoning. The GA, 

(v&oIPgOPgI^o) (^o|Q|^o) .29^ 

{^o\PgPg\^o) ^ («'o|*o) ' ^ ' 

approximates the expectation value within the projected state Pg |vI'o) by a 
corresponding statistical weight go multiplying the matrix element within the 
unprojected wave function |^'o). To determine the Gutzwiller renormalisation 
factor go we can either invoke Hilbert space counting arguments [20l 11201 1122] , 
or consider the limit of infinite dimensions {d = oo) [1231 11241 11251 In the 

following, we review both techniques and compare the respective results. 

^^To avoid confusion, in this section we denote density operators with a 'hat', and write, 
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3.1.1 Gutzwiller renormalisation factors by counting arguments 



Hilbert space counting arguments enable us to derive the renormalisation factor 
go through simple physical reasoning. We may use, 

9o « -r— (30) 

with 1^) = PgI^o), as defining the factor go\ (■•■)* denotes the expectation 
value with respect to the (projected) wave function j^f). In the GA, the ratio 
in ([50]) is determined by neglecting correlations in the wave functions |^) and 
|\l/o). The physical quantity which determines the theory is the occupancy at 
any site i. Thus, one calculates the probabilities for a site i to be empty, singly 
occupied with spin a, and doubly occupied, respectively. These probabilities 
are obtained by considering the Hilbert space restrictions and are summarised 
for 1^1') and |>I'o) in table [T] In this context, we should note that the densities 
before projection (n^, and n°|) and after projection (n^, riii and rii^) may 
differ. This is due to the projection operator, Pq = 11^(1 ~ "'il'^iT)j which 
can, e.g. remove more terms with an "f-clectron than a J,-electron on site i. 
Such effects become of importance for Gutzwiller projection in antiferromag- 
netic, charge ordered, or grand-canonical states. Keeping this caveat in mind, 
the expectation values in (j30p can be calculated approximately by considering 
the probability amplitudes of 'bra'- and 'ket'-configurations that contribute. We 
obtain the Gutzwiller renormalisation factor by calculating the ratio between 
these approximate expectation values. Although, wc neglect any off-site corre- 
lations in the derivation of the Gutzwiller renormalisation factor, the GA itself 
incorporates additional correlations by the expectation value of O in |^'o)- 
Extensions of the GA, which incorporate more correlation effects, were proposed 
by Ogata and Himeda ( |107| . see also section [4.4. 1[) and Hsu |152| . 

To illustrate the above scheme, we consider the expectation value of the 
hopping element, (cl^Cji). For a projected state, \^) = Pg\^o), we can write 

(cItCjt)* ^ ((1 - "a)clt(l - "ji)cjT)* ■ (31) 
We then perform the GA for the right hand side of (j3ip . which is written in 



occupancy on site i 


probabilities 


in *) 


in |*o) 


((l-na)(l-^«T)> 


1 - Jii 


(l-n^,)(l-n^,) 


(na(l - rijt)) 


nil 


n^i(l-n^T) 


("n(i - nti)) 


nil 


<T(l-<i) 







„o' 

^il^il 



Table 1: Probability for different occupancies on site i in j'l') and |^'o)- We 
distinguish between the densities before projection (n?, ti^^ and n°|) and after 
projection [ui, and n.i|). 
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Figure 18: Required bar- and ket-configurations, so that {{l—fiii)cl^ {l—hji)cj'\) 
contributes in |^) and |^'o)- Configurations that do not contribute to ((1 — 
~ S'l'e crossed out. 



It is convenient to 



terms of projected operators (1 — hi^)cl^ and (1 
rewrite the matrix elements in this manner before performing the GA, since it 
guarantees agreement with the infinite dimensions approach. Next we consider 
the probabihty for ((1 — ?^ij,)cj|(l — nj|)cj|) in j^I') and |^o)- Configurations can 
only contribute if the bra-vector has a single |-electron on site i and a vacancy 
on site j. For the ket-vector the interchanged occupancies are necessary, i.e. a 
single t-electron on site j, and a vacancy on site i. The corresponding hopping 
process is illustrated in figure[THl With the help of table[l]we find the amplitudes 
of the bra- and ket-contribution, and the product gives the probability in |^), 



and in |^'o), 



Kt(1 - rij)]^^^ ■ Kt(1 



11/2 



1/2 



(32) 



(33) 



The square roots stem from the fact that both bra- and ket-vector only provide 
amplitudes; the probability is obtained by a product of two amplitudes. 
Combining ^ and ^ yields 



((1 - na)cl|(i - ^ji)cjT)|*> 



((1 



^a)cl|(i 

9t 



1 



[n,|(l - nj)nj^{l - ?i»)]^^^ 

«°t(i-"°tKt(i-«°t)''^' 



(34) 
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The expectation value in j'J') is now obtained by renormalising the unprojected 
value by pi]) . 

(cItCjt)* = {{'^-riii)cl^{l-nji)cj^)^ (35a) 
« gt {il~n,i)cl^^{l-nji)cj^)^„ (35b) 
« 5t (cl|Cj|)*o . (35c) 



In the last row of ([55)) wc decoupled the densities in |^'o)- The Gutzwiller 
renormalisation factor is then. 



9t 



(36) 



We emphasise that the decoupling in p5c[) is controlled in the limit of infinite 
dimensions, viz. all neglected decouplings yield off-site correlations of higher 
ordei0 and thus vanish for d = oo. Violating this rule causes deviations from 
the mathematical thoroughness of the infinite dimension scheme. 

For the full determination of the Gutzwiller renormalisation factor in p6p . it 
is necessary to evaluate the dependence of the densities after projection relative 
to the densities prior to projection. The situation is particularly simple for a 
homogeneous wave functions with fixed particle number and spin symmetry, 
where n^^ = n^^ = = n/2 on each site i. Then, Ui-^ = ni[ = n/2, and the 
Gutzwiller factor simplifies to the well-known result, 

1—71 , ^ 

which incorporates the fact that the kinetic energy in j^*) is connected to the 
motion of holes, vanishing in the undoped case. 

However, the relation of the riig- with respect to the n'^^ become more subtle, 
if we consider, e.g. an antiferromagnet with sublattice magnetisation m, where 
'"■Act ~ "-/^ i m, and, n'\^ ~ ^b-ct (sublattices A and B, a In this 

case, riio- 7^ n^^, and we must invoke counting arguments to determine n.i^. We 
consider a canonical ensemble, where the overall particle density is the same 
before and after projection (rii = n^i +nii ~ H-tt-^^ — "-i* = n). Furthermore, 
the density riia is necessarily related to the probability of finding a single a- 
electron at site i in |5'o)- Thus, ni„ (x n°^(l — n^_^). Due the conserved particle 
density, 

JT-iT + ^il — n — , (38) 



Strictly speaking, we violate this rule by neglecting decouplings which include on-site 
pairing, (ct^c|^). However, we work in the fully projected limit, i.e. \'^) does not allow for 
on-site pairing. It is thus reasonable to prohibit on-site pairing in I'J'o) as well and to set 
(ct c+,> = 0. 
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and so, 

n,, = nt{l - ^t.)—£oZir ■ (39) 

Inserting tliis expression in tire numerator of p6p gives tlie Gutzwiller renor- 
malisation factor, 

^* = 1 - 2nOnO/n ' ^^^^ 

where and are from the same site. We note tliat ([5^ is valid for sites i 
and j on the same as well as on different sublattices as one can show easily and 
reduces to, gt = (1 — ~ '^/2), in the non-magnetic limit, n'^ = n/2. 

The situation becomes yet more complicated if we consider states with an 
inhomogeneous particle density, where it is difficult to determine n^, rui , and n^i . 
Therefore, most authors assume rii = n^. However, this assumption is incorrect, 
because the operator Pq ~ X]i(l^^jT ^-Uows for changes in the local particle 
density. An elegant solution is to redefine the operator Pq, so that Ui = or 
even ni„ = n^^. This conservation of local particle densities can be achieved 
by incorporating appropriate fugacity factors (that describe the local chemical 
potential) into a new operator Pq (Gutzwiller correlator), which is then not a 
projection operator any more. The redefined operator Pq still allows to present 
any projected wave function as j'l') — Pc\^o), however, the unprojected wave 
function |^o) wiU generally differ from I^I/o) defined by |*) = Pg|*o)- The 
use of Pq instead of Pq is often not explicitly stated in literature, although 
the assumed conservation of densities is only valid for a generalised Gutzwiller 
correlator Pq- Such a clear distinction between Pq and Pq becomes particularly 
important when results from the GA are compared to VMC calculations that 
implement the original Gutzwiller projector Pq. The non-conserving of local 
particle densities by the operator Pq also explains discrepancies between VMC 
calculations in the canonical and the grand canonical scheme [5]. This will be 
discussed in more detail in section [321 

Before turning to the d = oo scheme, let us discuss the the Gutzwiller 
renormalisation factor gs for the superexchange interaction, defined by 

(SjSj)* = .g5(SjSj)*o . (41) 



We first consider the GA for the contribution {S^ S 



I.e. 



(S+Sr)^ = (S+S-)^, . (42) 

The procedure resembles the derivation of gt. We note that the process S^S~ 

requires, an |-spin on site i and a |-spin on site j in the bra-vector, and the 
reverse in the ket-vector. Therefore, the probability becomes 

{n.^-njinnnj^y^^ (43) 

in the state j^E'), while it is 

[^(1 - n^xX(l - «?tX(1 - 4)"?t(1 - ]'^' (44) 
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in the state |4'o)- Using rii^a from ([M)) yields, 

± _ 1 

~ (1 - 2n<^7iO/7i)2 ■ ^^^> 

One can show again, that above formula results also for the case of sites belong- 
ing to the same sublattice. 

Next we evaluate the GA for the diagonal contribution to the superexchange, 

{S!S;h = 9s {S-S^h„ . (46) 

Here, we use, = l/2(ni| — fiii), and write, 

4 {S-S-) = (ni|(l - n^l)nj'^{l - n^j,)) + {n^l{l - n,t)nji(l - nj-^)) 

- (^*t(i - - "jt)) - - "jt)"jt(1 - "ji)> > (47) 

which is valid for any wave function. The Gutzwillcr approximations of the 
terms in (|47p give a common renormalisation factor, 

= (1 - 2n0n;/n)2 ' ^^^^ 

This is seen by considering the term — nix)nj|(l — fiji)) in (|47p . as an 

example. By applying the probabilities from table [U we obtain, 

,0 \„o n — nrr = 9s > l^yj 



where using ([55]) for Ui^riji directly confirms (pS]) . Since all density terms of 
(|47p renormalise in exact the same manner, gives the correct renormalisation 
factor for (S^S^)^ in 

From (|45p and (j48p . we find a common Gutzwiller renormalisation factor, 
5s = .9s — .9s : for (|41P - which simplifies to, 

^ (1-1/2)2 ' (50) 

in the non-magnetic limit, nj^ = n/2. At half-filling, n — 1 and gs 4, the 
magnetic correlations arc four times as pronounced in \^!) than in |^'o)- We 
note that Gutzwillcr approximations for other quantities arc easily obtained by 
following the same reasoning as for gt and gs- 



3.1.2 Gutzwiller renormalisation factors in infinite dimensions 

The effects of the Gutzwiller correlator can be evaluated exactly in the limit 
of infinite dimensions [1231 1124] . Gebhard [124] showed that a simple diagram- 
matic evaluation is possible for d = oo. Using the Gutzwiller renormalisation 
factors from d = oo for finite dimensions corresponds to a mean field approx- 
imation. Thus the d = oo approach provides a systematic way to calculate 
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Gutzwiller factors. Typically, one is interested in the doping dependence of 
such factors and here, the results from d = oo are in qualitative agreement with 
Gutzwiller factors calculated using counting arguments. Discrepancies between 
the two methods are merely quantitative. Here, we summarise the calculation 
of Gutzwiller factors in the limit d = oo, for fully projected states. The reader 
is referred to recent works of Biinemann, et al. [12511126] for a detailed account. 

To simplify calculations, the Gutzwiller projector Pq is reformulated as a 
Gutzwiller correlator Pq within the d = oo scheme. This redefinition agrees 
with the one discussed earlier and ensures that local densities are conserved, 
viz. Tiicr = n^^. The Gutzwiller correlator, Pq = Ili-PG.ii written as a 
product of local correlators, 

pG^^ = (1 - n,^){l - h,^) + \j n,^il ~h,i) + n.^il - n,^) . (51) 

Physically, the parameters A° and Af allow to weight locally the probabilities 
to find empty sites and sites occupied with a spin a, respectively. The A?, Xj, 
and Aj are determined by the constraints, 

{PD^o ^ 1, (52) 

{PG,^n^aPG.^)^,, = {fl^a) q, ^ = n°^„ . (53) 

Equation guarantees the normalisation, (*|*) = {^q\PgPg\'^o) = 1- of 
the projected wave function and equation (j53p provides the conservation of 
local densities. Evaluating these equations, we find. 



- V(l-na)(l-n,T) ' ^^^^ 
K = JjT^ ^- (55) 

y (1 - ni-a) 

Using these parameters in the Gutzwiller correlator Pq guarantees via ([5^ a 
conserved norm and via (|53p conserved spin densities for any projected wave 
function, l^*) = PqI'^'q). The GA for an operator Oij acting on the sites i and j 
is now obtained by neglecting all correlations except those between sites i and 
j. This procedure becomes exact in infinite dimensions and is written as. 



(*0|^^GOy^G|*0) = {^o\PGaPG.JO,,PG^^PG.J\'<^^o) ■ (56) 

Decoupling the right hand site and neglecting all off-site correlations of higher 
order, provides the exact solution for d — oo, which agrees with the results from 
counting arguments presented in section [3. 1.1 1 

As an example, we consider the hopping process, (c|-|.Cj|) p^i^^^- Using (|56[) . 



41 



wc find, 



(*o|Pgc1^c,tPg|*o) = {^o\PgA^Pg,,Pg,jC,^Pg.j\^o) (57a) 
= A|aOaOa] (^-oKl - n,x)ct^(l - n,i)c,T|*o) (57b) 
AjA°A°Aj(l - n,J(l - (*o|c,V,T l^-o) , (57c) 



where we decoupled the densities in the last row as already done in the discussion 
using counting arguments. Equation (j57p is exact in infinite dimensions, and 
gives the Gutzwiller renormalisation factor, 



(1 - nj){l - n^) 



V (l-'^iT)(l-"n) 

which agrees with (|36p if we assume locally conserved densities. However, we 
note that this result differs from which incorporates the changed spin 

densities due to the projection operator Pq. The scheme presented above is 
applicable to any kind of operator, and gives the exact result for d — oo. It 
provides an useful check for results derived from counting arguments. Never- 
theless we must keep in mind that results may differ depending on our choice 
of Pg and Pd'i'o) or Pq and Pg\'^o), in the counting arguments leading to the 
derivation of the Gutzwiller renormalisation factors. 



3.2 Gutzwiller approximation in the canonical and the 
grand canonical scheme 

In this subsection we follow Edegger, et al. [5] and study the effects of projection 
on superconducting BCS wave functions, 

|vI'o) = |BCS)sn("k + t'k4TC^-ki)|0) . (59) 

k 

Since |^'o) = |BCS) exhibits particle number fluctuations, the projection oper- 
ator Pq can change the average particle number N of the wave function, i.e. in 
general, 

(vl/ol^I/o) ^ (^-olPai^o) ■ ^ ^ 

In the above equation, the equality between the l.h.s. and r.h.s. could be 
recovered by replacing the Gutzwiller projector Pq by a Gutzwiller correlator Pq 
which conserves local densities as discussed in the previous subsection. Here, we 
follow a different route to compensate the effects of projection by using a fugacity 
factor in the wave function. This ansatz explains differences observed between 
VMC calculations in the canonical framework (fixed particle number) and the 
grand canonical ensemble (fluctuating particle number) using the corresponding 
GA. 
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3.2.1 Incorporation of a fugacity factor 

We first examine the particle number distributions to illustrate the effect of 
the projection operator Pq in the projected Hilbert space. Towards this end, 
we write the average numbers, N^^") (N) in the unprojected (projected) Hilbert 
space, as, 



7V(")=^iVp(^', TV^^iVp^, (61) 

TV N 



where. 



(0) _ (^olPjvl^o) _ {^o\Pg Pn Pg I^o) 

" (*o|*o) ' {^o\Pg Pg I^Po) ' ^ ^ 

are the particle number distributions in the unprojected and projected BCS 
wave functions respectively. Here, the operator P^r describes the projection 
onto terms with particle number N. As discussed in [5], we can relate the 
particle number distributions before and after projection by 

(^-ol Pg Pn Pg |*o) (*o| Pn |*o) 

9N UT, X , (63) 



(^-oIPg Pg I^-o) (^-01*0 



PN 



with 

^ (^ol^o) (^o|Pg Pn i^cl^o) . .. 

(v1/o|PgPg|*o) (^oIPatI^o) ■ ^ ' 



—G(— const) 

The above equation ([551) describes the GA for the projection operator P/v with 
the corresponding rcnormalisation factor gi^. The parameter C is an irrelevant 
constant, which does not depend on N . It follows that if we were to impose 
the condition that the average particle numbers before and after projection be 
identical, a factor g~^^ needs to be included in (j6ip . Then, from (|6ip and ([63]) . 
we obtain the particle number after projection Nnow, 



-I (0) 

7V_^ y^N^p^ = ViV^^ = 7V(«) , (65) 

N 3N N 3N 

which is the desired result. 

This procedure can be implemented for the wave function l^fo). Since the 
BCS wave function is a linear superposition of states with particle numbers 
. . . , — 2, iV, TV + 2, . . ., we consider the effect of projection on two states whose 
particle numbers differ by two. Then, the ratio is 

.2 _ 9N+2 f L- N Y f '^^'^ V 2 

where the factors gN were evaluated combinatorially in the thermodynamic 
limit [5]. We note that the fugacity factor / is equal to gt, the Gutzwillcr 
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renormalisation factor for the hopping term (see also section [57T|) . Equation 
((66|) shows that the projection operator acts unequally on the N and N + 2 
particle states; the renormalisation of the weight of the + 2 particle states 
gN+2, is gf times the weight of the N particle states, g^- This effect can be 
rectified as in ([65]) by multiplying every Cooper pair cjj.|cLjj.| by a amplitude 
in the BCS wave function. Alternatively (following Anderson [153] ). we can 
multiply every empty state by the factor gt, and write 

\^i'')=U ^ 7 10). (67) 

Then again by construction, the fugacity factor f = gt in ([57]) ensures that 
the projected wave function Pg\'^q ') and the unprojected wave function \^o) 
without fugacity factor have the same particle number. The denominator in 
(j67[) is the new normalisation factor. Two points deserve further attention. The 
first is that a relative phase factor between N and N + 2 particle states in the 
projected BCS wave function can be absorbed into the definition of the fugacity 
factor. The second point is that our ansatz for the fugacity factor assumes that 
the effects of the projection operator are independent of fc; i.e., the fugacity 
factor we obtain is independent of k. We do not sec a priori why the fugacity 
factor cannot depend on k. This is an interesting line of investigation since a k- 
dependent fugacity factor in the Gutzwiller - BCS wave function would lead to 
(experimentally) verifiable consequence. 

3.2.2 Singular particle number renormalisation close to half-filling 

We saw that the inclusion of the fugacity factor is necessary for the average 
particle number in a BCS wave function to remain unchanged when projecting 
out all doubly occupied sites. Alternatively, one might ask what is the effect of 
the projection operator on a BCS wave function, if projection changes the mean 
particle number, because no fugacity factor is introduced. 

In this situation, as shown by Edegger, et at [5], the the particle density 
after projection is determined by the self-consistent equation. 



which can be solved iteratively with gt specified by the particle number in A^after- 
Since the particle density in the state \'^o) before projection is given by, 

— A'bcforc 2 ■r-^ 2 /•/;n\ 
"^bcforo = = ]^ l^kl , (69) 

k 

equation 168]) provides a way to calculate the particle number in the state j'l'o) 
after projection, whenever the particle number in the state j'l'o) before projec- 
tion is known as a function of the identical factors and v^. 
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Figure 19: Particle density before projection nbcforc, equation (|69p . as a function 
of the particle density after projection riafter, equation (j68p . for different d- 
wave order parameters A. The dashed line indicates the Fermi liquid result 

"^bcforc = J^aftor- FrOm [5]. 



In the following, we discuss numerical solutions of (j68p . where we use the 
standard BCS expressions for a d-wave superconductor, 



2 

ft = 



with, 



i?k = VA^+e^, (71a) 
Ak = A ( cos fcj; — cos fcj, ) , (71b) 
= —2 {cos kx + cos ky) — ^ . (71c) 

Therefore, the only free parameters, which must be specified within the calcu- 
lations, are the chemical potential /i and the order parameter A. 

The particle numbers (before and after projection) for fixed values of the 
order parameter A can now be determined as a function of the chemical poten- 
tial. The results for the particle densities are shown in figure [T21 The results 
clearly show that the particle density before projection attains its maximal value 
(jT-bcforo = 2), when Tiafter = 1 (half- filling) . This result holds for any finite value 
of the order parameter A. The case of half-filling, riafter — *■ 1, is therefore sin- 
gular in the grand canonical scheme and large deviations with respect to the 
canonical framework can be expected. In the opposite limit, viz. low densities 
of electrons, nboforo converges to the value of riaftcr as expected. As illustrated 
by the results in figure [191 the size of the intermediate region depends on the 
magnitude of the order parameter A, i.e. the effects of projection are largest 
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Figure 20: The particle density after projection riaftcr as a function of the pa- 
rameter A for a d-wave BCS state at various chemical potentials /i. The figure 
shows a comparison between results from equation (|68p (solid lines) and the 
VMC results of Yokoyama and Shiba [57] [for 6x6- (circles) and 8x8 -lattices 
(squares)]. Numbers in the figure denote the chemical potentials of the corre- 
sponding curves. From [5]. 

for wavcfunction with a large A and consequently large particle number fluctu- 
ations. 

To check the accuracy of equation ([55)1 we compare it with VMC calculations 
of Yokoyama and Shiba (YS) [97j, who numerically studied a projected BCS 
wave functions with fluctuating particle number (but without a fugacity factor). 
YS determined the particle density of the projected d-wave state Pq I^'o) as a 
function of the parameter A for various flxed chemical potentials /i within a 
grand canonical scheme (see VMC data in figure EO)) . Since the unprojected 
wave function was specified through ([70|) - (|71|) in the VMC calculation, we 
can also determine the relation between riaftcr and A by (j68p . As shown in figure 
[20l the results from ([68]) are in good qualitative agreement with the VMC data 
of YS. Small discrepancies arc mostly explained by finite size corrections in the 
VMC calculation (VMC calculation only for 6 x 6 and 8 x 8-lattices). Figure 
[20] clearly reveals the singular effect of the projection near the insulating phase 
(half filling), where the chemical potential diverges to infinity. 

3.2.3 Gutzwiller renormalisation factors in the canonical and the 
grand canonical ensemble 

Next wc discuss the differences between the Gutzwiller approximation in the 
canonical and grand canonical scheme. The validity of the analytical expressions 
derived in this section can be confirmed [5] by a comparison with numerically 
exact VMC calculations [T71I97]. 
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Figure 21: (a) The kinetic energy i?'^' and (b) the energy of the remaining 
terms i?^^-* per site of the t-J model as a function of A for the d-wave state 
at a fiUing n = 0.9. Fixed particle (can., circles) VMC data [T7] and grand 
canonical (grand can., squares) VMC data [97] are compared. The dashed/solid 
lines represent the corresponding Gutzwiller approximations (GA). From [S]. 



We first consider the canonical case, where we are interested in the ex- 
pectation value of an operator O calculated within a projected wave function 
PnPg \ ^'o) with fixed particle number. The GA corresponding to the operator 
O is given by, 



{^o\Pg Pn O Pn PgI^-o 



(*o|Pg Pn Pal'f 







(72a) 



with the Gutzwiller renormalisation factor go and the projector P/v onto terms 
with particle number A'^. The term (j72ap represents a quantity which can be 
calculated exactly by the VMC scheme with fixed particle number [T71I2I] • Since 
the particle number is fixed, the usual Gutzwiller approximation can be invoked, 
leading to (j72bp . The equality to the last row is guaranteed only when A^ is 
equal to the average particle number of | ^^q) [N = N) . Under this condition a 
transformation from a canonical to a grand canonical ensemble is valid in the 
pre-projected Hilbert-space. 

In the grand canonical scheme the expectation value of O is calculated with 
a particle number non-conserving wave function. Therefore, this scheme must 
be modified as follows, 

{^If^lPaO Pcl^lf) _ (MOl'^o) 



where Pg \ ^o''^'') is the projected c?-wave state including a fugacity factor ([S7|) 
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Figure 22: Comparison of the canonical and grand canonical scheme at filling 
(dashed lines) n = 0.9 and (sohd lines) n = 0.99. (a) The kinetic energy 
and (b) the energy of the remaining terms E'-^'> per site of the t-J model are 
shown as a function of A for the d-wave state. The results are obtained by 
Gutzwiller approximations, i.e. via (|72p and (j73p for the canonical and the 
grand canonical scheme, respectively. The calculations for the grand canonical 
scheme follow the steps given by [5]- 



This correction is essential to guarantee the validity of the Gutzwiller approx- 
imation, because without it, the left hand side (Ihs) and the right hand side 
(rhs) of ([73]) would correspond to states with different mean particle numbers. 
Equation ([73]) shows that a fugacity factor must be included into wave functions 
with particle number fluctuations when used for a Gutzwiller approximation. 
A comparison of ([7^ and (|73p gives the important result. 

This equation shows that to obtain identical results in the grand canonical (rhs) 
and canonical (Ihs) scheme, one has to use different wave functions. The wave 
function | ^'g ) is a d-wave state including a fugacity factor, whereas | ^^/q) is a 
pure d-wave state. 

The arguments leading up to (|73p and (j74p can be verified by a compari- 
son with VMC studies. In figure 1211 we show Gutzwiller approximations from 
Edcgger, et al. [5] together with the corresponding VMC calculations, i.e. fixed 
particle number VMC [17j for the canonical scheme and particle number non- 
conserving VMC [97] for the grand canonical scheme. The figure shows that 
canonical and grand canonical approaches yield different energies, however the 
GA qualitatively matches the corresponding VMC results. The differences in 
the two schemes are due the projection operator Pg, which changes the particle 
number in a grand canonical wave function. For these two methods to yield 
the same results, a fugacity factor must be incorporated in the wave function 
when working in a grand canonical ensemble. We note that the discrepancies 
between the canonical and the grand canonical scheme increase significantly 
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towards half-filling as illustrated in figure This effect is due to the strong 
renormalisation of the particle density in this limit (see figure [19]). 

3.3 Gutzwiller approximation for partially projected states 

Thus far, we discussed the Gutzwiller renormalisation scheme for fully projected 
wave functions. It is however necessary to consider sometimes, partially pro- 
jected states of the form, 

- Pi 1*0), Pi =l[{l- n^^n^i) ■ (75) 

The wave function describes a state where double occupancies are projected 
out on all sites except the site I, which we call the "reservoir" site. The reason 
for the appearance of reservoir sites can be seen as follows. Consider, for exam- 
ple, the operator PgCi-\, which can be rewritten as c^Pl- Such commutations 
become necessary, e.g. for the calculation of the quasiparticle weight (discussed 
in section [6]), where partially projected states arise inevitably. 

Before discussing \^[) in more detail, we remark that the notation 'partially 
projected' is also used for a projection operator, 

Pa = ]^(1 -afii^hii) , (76) 

i 

with a S [0, 1]. The operator Pa is used for studying Hubbard-like model with 
'partially projected' wave functions PaY^o) (see also section l3J|) . Here, the 
parameter a controls the total number of double occupancies, whereas P/ in 
(|75)) yields a fully projected state with only a single unprojected reservoir site 
i. We emphasise that the respective Gutzwiller approximations for these two 
projection operators are fundamentally different. 

Below we follow the work of Fukushima, et al. [5], who developed an ana- 
lytical method to calculate expectation values for partially projected states [as 
defined in ([751) ]. The calculations rest on counting arguments, however, similar 
results can in principle be obtained within the infinite dimensions approach. 
We first determine the local occupancy of the reservoir site, which is then used 
to derive the Gutzwiller renormalisation factors of specific expectation values. 
We also provide a comparison to VMC calculations to test the validity of the 
approximation and determine its limitations. 

We are interested in expectation values such as 

mm) (*oi«'o) ' ^ ' 

that generalise the Gutzwiller approximation to partially projected wave func- 
tions. Note that the reservoir site does not have a special role in the unprojected 
wave function I'I'o)- This is in contrast to the impurity problem, where an im- 
purity site would break the translational invariance of both the unprojected and 
of the projected wave function. 
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3.3.1 Occupancy of the reservoir site 

The GutzwiUer approximation in ([77)) can be performed by counting arguments 
as in the fully projected case (see section [?.l.ip . However, the occupancy of the 
reservoir will differ from the occupancy of a fully projected site, an effect that 
must be considered when deriving GutzwiUer renormalisation factors. 

Fukushima, et al. [6] showed that the probabilities for the reservoir site to 
be empty, single occupied, or double occupied are, 



((1 - fi;|)(l - hii))^ 



respectively. Here 



= Xricr ? 
= 1-X 



(1 



(l-nT)(l-nx) 
(1 - n)ncr 



(l-^T)(l-"i) ' 



X 



(^oIPg^gI^o) 
(vPolP/P/lvI/o) 



(78) 

(79) 
(80) 

(81) 



is defined as the ratio between the normalisations and This ratio 

can be estimated by a GutzwiUer approximation [S] and becomes 

1-71 



X 



(l-nT)(l-ni) 



(82) 



in the thermodynamic limit, where n„ = N^/L {a = t, i) and n = + are 
the respective particle densities. We note that X vanishes at half-filling. Conse- 
quently the reservoir becomes exactly doubly occupied, i.e. lim„_,i — 1. 

3.3.2 Renormalisation of mixed hopping terms 



The occupancies of the reservoir site, (|78 l) -(|80 l) . directly enter the respective 
GutzwiUer renormalisation factor gQ. We consider here, as an example, the 
mixed hopping term. 



(*ol^/cLc™.P/|*o> ^ , (*o|cLc^.l*o) 



(*o|P/P/|*o 



9t 



(*o|^'o) 



(83) 



where I denotes the reservoir site and m ^ I is a fully projected site. Following 
the arguments leading to (j3ip . we rewrite. 



(c||Cm-f)3(J — (c|| (1 — 7T.,„|)c,„|),f J 

~ 9't (c||(l - fimi)CmT)*(, 



(84a) 
(84b) 
(84c) 



=9t 
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Figure 23: Required bra- and ket-configurations. so that {c^(l — hjni)cm'i) 
contribute in when I is a reservoir site (indicated by a shaded background). 
The last row presents the sum from the two possible contributions as used in 
([85|) . Boxes with white background indicate the fully projected site m. 



As for ([3T|) . we perform the GA for the rhs of (|84ap to guarantee agreement 
with the infinite dimensions approach. However, the decoupling of (1 — Umi) in 
(|84c|) becomes exact in infinite dimensions and we can recover the GA of ([83|) . 

In analogy with the calculations in section 13.11 we consider the probability 
for (4^(1 - njni)cmi) in l^'J) to determine the corresponding Gutzwiller renor- 
malisation factor g[ entering in (|84cp . We illustrate the two configurations that 
can contribute, together with the resulting probability from combining the bra- 
and ket-vectors, in figure [23l Using ((78|) - ([80|) for the partially projected site 
(grey in figure 1^51) . we find the probability. 



(^[Xni^X{l - ni)f^ + [(1 - • [(1 - nm)nm^\ 



1/2 



, -1 l/2\ 

= X{[ni^il-ni)f'+ ^^ni^l \ ■ [{1 - nr„)n^^f\ (85) 



in l^'J). With njo- = ~ n„ and (|82p . the above expression simplifies to, 

X [n,{l - n) + nl] ^ X n^{l - n,) . (86) 
For the respective probabilities in |^o)i we use table [T] and obtain, 

Kt(1-0]'^'- KT(l-"™i)(l-'Ci)(l-<.T)]'^' . (87) 

As pointed out in section 13.11 Ji^^ = n^n^ ~ no-, for non- magnetic wave 
function. We then get the renormalisation factor g'^ from the ratio of (|86[) and 
(EH), i.e. 
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hole doping 



Figure 24: Double oceupancy of the reservoir site (c?)^j = ('^iT^-u)*; a func- 
tion of doping, for the partially projected Fermi sea. Note the good agreement 
between the Gutzwiller result (solid line), equation ([80|) . and the VMC results 
for the projected Fermi sea (open circles). Statistical errors and finite-size cor- 
rections are estimated to be smaller than the symbols. From [5]. 



Together with (|84cp . we obtain the renormalisation factor, 

1 - n 

g[^{l- n^yg't ^X = ^ , (89) 

for the GA in ([551) . Other expectation values in partially projected states can 
be calculated similarly. In section [6l we will use the same scheme in calculating 
the quasiparticle weight. 



3.3.3 Comparison of the GA for partially projected states with VMC 
calculations 

Before concluding this section, we illustrate a comparison between ((50)) and 
VMC results for = (^/T^u)*; = 1 — V. Fukushima. et aL[6] found that 

the results obtained by a generalised Gutzwiller approximation are in excellent 
qualitative agreement with the VMC results for a partially projected Fermi 
sea as shown in figure [24l We also used VMC to obtain the same quantity 
using projected s/d-wave BCS stateJ^ as variational states in the simulation. 
The results for {d)^^'^ in BCS states are shown in figure [25l In contrast to 
the projected Fermi sea, a clear deviation from the Gutzwiller approximation is 
seen. This underscores the importance of pairing correlations in the unprojected 
wave function that are not completely taken into account by the Gutzwiller 

"The BCS states are defined by, |t)kp = 1/2(1 -5k/£^k), and UkU* = Ak/(2£;k), where 
^k = — 2 (cos fca: -I- COS ky) — /i and i?k = ^|AkP + ^k [s-wave: Ak = A, d-wavc: Ak = 
A (cos kx — cos ky)]. 
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hole doping 

Figure 25: Double occupancy of the reservoir site (c?)^j = ('^iT^-u)*; ^ func- 
tion of doping, for the partially projected BCS wave function. The solid line is 
the GA result from ([SO]) . The parameterisation follows [147j . Statistical errors 
and finite-size corrections for the VMC results are estimated to be smaller than 
the symbols. From 0. 

approximation scheme. It also explains to a certain extent, the discrepancies 
between the VMC calculations and the GA for the quasiparticle weight as we 
will discuss in section [T] 

To clarify the limitations of the Gutzwiller approximation for projected su- 
perconducting states in more detail, VMC can be used to calculate the hole 
density in the vicinity of the reservoir site. In the half-filled limit the reser- 
voir site is double occupied and therefore a single hole is distributed among the 
remaining fully projected sites. The VMC calculations of Fukushima et al. 
[5] show very different density oscillations for the projected Fermi sea and the 
projected d-wave state. 

Figure UHl shows VMC results for the hole density 

nh{m) = (1 - , 

in the partially projected state The sites m are distinct from the reservoir 
site I (marked by a cross in the figure). All results shown correspond to half 
filling; viz. n-f = = 0.5. We choose A = 1 for the projected BCS d-wave 
state. For the Fermi sea, we see that the hole is distributed more uniformly than 
in the d-wave case even though the diagonal direction has a larger probability of 
being occupied by a hole. The d-wave has a quasi checker-board pattern where 
only one of four sites is black, and the hole tends to be near the reservoir site. 
The VMC results for the projected BCS wave functions are strikingly different 
in that the hole density is not uniform. On the other hand, the Gutzwiller ap- 
proximation would be exact, if all states in the Hilbert space contribute equally 
to the wave function. That would correspond to a uniform density of holes. 
Clearly, some limitations in the Gutzwiller approximation show up when treat- 
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5 10 15 20 25 30 35 5 10 15 20 25 30 35 

Figure 26: VMC results for the hole density nh{m) = ((1 — rim))^'^ (colour 
coding: white/black correspond to high/low values of nh{m)) in the partially 
projected state for sites m other than the reservoir site / (marked by the 
cross). Left: Fermi sea. Right: d-wave state. From [6] 

ing projected superconducting wave functions. This is in agreement with our 
previous considerations, where we found that the functional form of X [ ([82|) . 
derived using Gutzwiller approximation] agrees well with the VMC calculations 
only for the projected Fermi sea, but not for the projected BCS state (see figure 
[Hand figure [13. 

4 Renormalised mean field theory: Basic ideas 
and recent extensions 

On the basis of the GA, Zhang, et al. [lOj derived a renormalised mean field 
theory (RMFT) for the t-J model. In this section, we present an overview on 
this approach, which plays a central role within Gutzwiller-RVB theory. We 
illustrate successes and recent extensions of the RMFT for the HTSC, derive 
the RMFT gap equation, and review its applications to the Hubbard model 
in the strong coupling limit. Further extensions to antiferromagnetic and in- 
homogeneous phases provide a quantitative description of the Cuprate phase 
diagram. 

4.1 Overview on the RMFT method 

We start with an overview and discuss how RMFT allows for a systematic 
treatment of the Hubbard Hamiltonian in the strong-coupling limit. We present 
the basic concepts in this subsection, and the method itself will be discussed in 
detail in the corresponding subsections. Figure [27l summarises the main steps 
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Figure 27: Schematic illustration of the RMFT method; see text for a detailed 
description. 



necessary for the strong-coupling treatment of the Hubbard model within the 
RMFT. In the following we refer to the individual steps illustrated in figure \T7\ 

As shown in figure 1271 the first step is to apply a canonical transformation 
e~*'^ to the Hubbard Hamiltonian removing hopping processes that change the 
number of doubly occupied sites. Doing so, we obtain the t-J Hamiltonian, 
which is defined in the subspace excluding double occupancy. The t-J Hamil- 
tonian provides an effective low energy Hamiltonian for the Hubbard model in 
the strong coupling limit as already discussed in detail in section 12.41 

Next we invoke the Gutzwillcr approximation to remove the restriction to 
projected states within the t-J Hamiltonian. As we will discuss in scction R:.2.H 
this procedure results in a renormalised Hamiltonian with terms weighted by 
the corresponding doping-dcpendcnt Gutzwillcr renormalisation factors (see also 
section [5T|) . 

We then perform a mean field decoupling for the renormalised Hamiltonian, 
focusing on hopping amplitudes £,,. = J2a {^lcr'^i+r<y)^o ^^"^ pairing amplitudes 
Ar = (c|^| c|^_|_j,| — c^|c|^_|_^|)>i,j, . Following this way, we find self-consistent gap 
equations for the mean field amplitudes (see section l4.2.2p . 

Solving the gap equations provides us with the mean field ground state 
l^'o) of the renormalised t-J Hamiltonian in the 'prc-projectcd' Hilbert space. 
Once l^o) is known, we can construct the Gutzwillcr projected state \'^) = 
Pg\^o): which then provides an approximate wave function for ground state of 
the projected t-J Hamiltonian. We can control the RMFT by using j^E') as a 
projected trial wave function within the VMC technique (sec section [S]). Thus, 
the above scheme provides a consistent framework to study Gutzwillcr projected 
wave functions by a combination of RMFT and VMC methods. 

The projected wave function j^*) allows for the calculation of relevant physi- 
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cal quantities as well as for the definition of excited states within the t-J model. 
To determine observables within the Hubbard Hamiltonian, we can employ the 
re-transformed wave function e~*'^|^') for the calculation of expectation values. 
Wc will discuss this approach in section 14.31 



4.2 Derivation of the RMFT gap equations 



In this subsection we review the work of Zhang, et al. [20j and develop a renor- 
malised mean field theory (RMFT) for the t-J model based on the Gutzwillcr 
renormalisation scheme (= GA, see section [O]) . To illustrate the RMFT, we 
start with the simplest form of the t-J Hamiltonian, 



t-J 



= Pg 



Pg 



(90) 



We restrict ourselves to nearest neighbour hopping t, and a superexchange in- 
teraction J. We neglect any further hopping parameters as well as additional 
contributions in the Hamiltonian like the density-density term and the correlated 
hopping terms, see equation ([6]). The effects of such extensions are discussed 
in section l473l where we consider an RMFT for the Hubbard model, including 
next nearest neighbour hopping matrix elements. 



4.2.1 Derivation of the renormalised t-J Hamiltonian 

Two steps are necessary to obtain explicit analytic expressions for the ground 
state of the t-J model (|90l) for various doping levels x, where x = 1 — n. The 
first is the Gutzwillcr approximation, where the effects of the projection Pq 
are taken into account by appropriate renormalisation factors. We search for a 
Gutzwillcr projected state Pcl'i'o) that minimises the energy expectation value, 



En = 



-t ^ ^ 



(*o|-Pg-Pg|«'o) 

By invoking a GA for (|9ip. we get rid of the projection operator Pq and obtain, 

The GA for the hopping term [first term in (jM]) ] has a renormalisation factor, 
gt = (1 — ^^)/(1 ^ '^/2), which was derived in the previous section, equation 
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(|57|) . For the superexchange term [second term in ([5^ ]. the renormalisation 
factor is = 1/(1 — •n/2y' , where we assume a homogeneous state without any 
sublattice magnetisation, see ([5n|. 

We may now determine the variational ground state by searching for the 
state l^o): that minimises the renormahsed t-J Hamihonian, Ht-j, defined as, 

Ht-j = -9tt (4,aCj.a + cl^c,.„) + 5s J ^ S, S, . (93) 

Once 1 4*0) is known, wc may consider the projected state, PgI^o): as a trial 
ground-state of Ht-j. 



4.2.2 Mean field decoupling of the renormalised Hamiltonian 

The next step in the derivation of the RMFT, see section HTTl is the realisation 
that Ht-j allows for several types of molecular-fields [71[5D]: For simplification 
we only concentrate on the singlet pairing amplitude, 

A,= {cl^4+^^~clcl^^^)^, , (94) 

and the hopping amplitude, 

1^ = ^(c|^c,+^<,)*o , (95) 

where r = x,y = {I, 0), (0, 1) connects nearest neighbour sites. This decoupling 
scheme of the renormalised Hamiltonian leads to a BCS ground state, 

l*o)=nK + ^^k4^cLk^) |0) , (96) 

k 

with, 

vl-Ul-^), (97) 



and "f^k ~ 1 ^ ^k- T^^'^ resulting gap equations are 



Ar = l/L^cos(kr)Ak/£;k , (98) 

k 

e; = -l/L^cos(kr)Ck/i^k , (99) 

k 

together with the condition, x ~ l/L^j^^k/^-k, for the hole-doping concentra- 
tion. The dispersion of the mean field excitations is given by, E'k = \/^k + Aj, 
where 

Ak = ( Ax cos + Ay cos ky) (100) 



^k = - ( 2gtt + cos k, - [2gtt + ^^iy ) cos ky - fi . (101) 
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doping X 

Figure 28: Doping dependence of the rf-wavc pairing amplitude A, the hopping 
amplitude ^, and the superconducting order parameter $, see (|102p . in the 
d-wave ground state for the t-J model ([M]) with J = The results were 

obtained self consistently by solving the RMFT gap equations. 

Equation ((98)) resembles the usual BCS gap equation, except that we consider 
independent pairing along the x- and the y-direction. Together with (j99p . we 
have four coupled gap equations (for A^;, Aj,, S.y), which must be solved self- 
consistently. The Ak is obviously related to pairing in the state |^'o), however, 
it is not identical to the superconducting order parameter in PqI'^o) as will be 
shown below. The becomes the renormalised dispersion in the absence of 
pairing and includes a chemical potential fi to regulate the particle density. 

4.2.3 Solutions of the RMFT gap equations 

The gap equation can be solved numerically. We present results obtained for 
J/t = 1/3, which is a reasonable choice for HTSC. However, we emphasise that 
the results presented below are not sensitive to this particular choice of J, i.e. 
the results stay quite similar for J/t = 0.2 — 0.5. We find that a c?-wave pairing 
state is stable for x < 0.35. In this case, A = jA^,! ~ |Aj,| with A.^ ~ ~Aj, and 
£^ = = . We illustrate the doping dependence of these quantities in figure 
[28l The superconducting order parameter, 

'^^\{4^4+rl-4A+r^)^\ > (102) 

is an expectation value in the projected ground states, j^*) = PqI^o), where r 
is a neighbouring site. Evaluating $ by the GA (section [33]) one finds that $ is 
renormalised as the hopping amplitude by gt, namely $ « gt A. As illustrated 
in figure 1281 $ vanishes linearly near a: = 0, while A continuously increases 
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towards half-filling. These results are in good agreement with VMC results 
|117[ 11471 \21\ and the experimentally observed Tc for the d-wave pairing in the 
HTSC. 

For the renormalised Hamiltonian, the i?k corresponds to the dispersion 
of the Bogoliubov quasiparticles, IvP^o) = Tktrl^o), with cr =ti j where the 
corresponding Bogoliubov operators are defined by, 7^^^ = '"kcLkx+'^kCkTj find, 
7j[.| = UkC^i — UkC_kj,, respectively. However, E'k also describes the excitation 
energy of the corresponding projected Gutzwillcr-Bogoliubov quasiparticles, 

|«'k.>^/^Gl*k,0> = ^G7Ll*o> . (103) 

To see why, one evaluates the expectation value of the t-J Hamiltonian with 
respect to l^ko-)- Because, l^'ko-) = -Pcl^'ko)' renormalised exactly as, 
l^f) = PgI^'o), we recover the renormalised Hamiltonian i/j-j, (|93p . by in- 
voking a GA. The state |^'k q) is now acting onto Ht^j, yielding in mean field 
decoupling a Bogoliubov quasiparticle with excitation energy i?k- Therefore, 
the gap Ak in ii'k = Vsk+^k corresponds to the quasiparticle gap in the 
projected superconducting state and is directly proportional to the mean field 
amplitude A in figure [251 We note that RMFT [5D] correctly reported the 
doping dependence of the d-wave gap, i.e. an increasing gap with decreasing 
doping, well before this behaviour was experimentally established. 

The above calculations follow the original work of Zhang et al. [20] Note that 
the results are restricted to a homogeneous and non magnetic phase. Therefore, 
the results cannot adequately describe the antiferromagnetic region of the phase 
diagram near half- filling as well as inhomogencous phases observed in HTSC. 
However, we emphasise that the Gutzwiller - RVB being a variational mean 
field theory can be extended to study such phases as well. In the following 
subsections, we describe some attempts made in this direction, that provide a 
more detailed description of the phase diagram of the HTSC. 

4.2.4 Local SU{2) symmetry in the half-filled limit 

In the limit of half filling, the kinetic energy renormalises to zero, since gt 
as X — > 0. The t-J model reduces to the antiferromagnetic Heisenberg model, 
which is conserved under local SU{2) gauge transformations pO l I154j . 

4t ^ "»'^It + -^'^'i ' ^^i- ~^ ~l^i4^ + ^i^n 1 (104) 

where a^a* -I- (3if3* = 1. The invariance of the Hamiltonian is due the spin 
operator Si, which is invariant under SU{2) transformations, as can be proved 
by applying (|104[) to the operators Sf and . For we find, 

St = clt^'i ^ ("''^It + Ac»i)(-/3*cJ^ + a*c,i) 

^ a,a* S+ + A/3* S+ = S+ . (105) 

The invariance of S~ and S'f under (|104p can be shown analogously. 
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Owing to the local 5/7(2) gauge symmetry, the renormalised mean field 
Hamiltonian has a large degeneracy in the representation of ground states at 
half-filling, as may be seen by transforming the mean field amplitudes and 
£,r, ((M|) and ([M)) . under (|104p . Some of the resulting (degenerate) states, that 
arc related to each other by SU{2) transformations, are summarised in table 
[21 Another example among the degenerate states in the SU (2) manifold is the 
staggered tt-Aux state [154] . 

eij = |eopcxp(i(-l)^^+^«J), (106) 

with a complex hopping amplitude = (cj^Cio-)o, but a vanishing pairing am- 
plitude Aij = (c^Ci|)o = 0. 

It is important to note that above degeneracy is not true in terms of the 
projected wave function, since it only results from using an undcr-dctcrmincd 
representation. In other words, the states that arc degenerate arc the unpro- 
jected states |^'o), but not the physical states PgI^'o) [20]. Therefore, the entire 
set of degenerate grounds states in the renormalised mean field Hamiltonian 
correspond (modulo a trivial phase factor) to a single projected state; for a 
proof see [20] . 

All states listed in the table[2|and p06p have the same superexchange energy, 
even at finite doping, due to the SU{2) invariance of this term. However, the 
kinetic energy T, ([8]), and the 3-site term H^^, (|10p . which only vanish at half- 
filling, are not invariant under the SU{2) transformation p04p . Therefore the 
degeneracy in j'l'o) is lifted at finite doping, where the d-wave pairing state is 
selected due to its lower kinetic energy. 

The SU{2) gauge symmetry of the superexchange term led to the speculation 
that at finite doping, when the degeneracy is lifted, some among the 'degenerate' 
states may compete with the d-wave state (also one of the degenerated states at 
half-filling). In particular, it was argued that staggered flux states could serve as 
a 'competing' and/or as a 'normal' state in the undcrdoped regime of the HTSC 
P^ni 1142] . Besides this competing order scenario, a spin-charge locking mecha- 
nism |113| resting on the presence of the SU(2) gauge symmetry at half-filling 
was suggested. The consequence of the degeneracy in the unprojected mean 



(i-wave pairing 


A^ = 


= -A, = 1, = 4 = C/V2 


d-wave density matrix* 


A^ = 


= Ky = = — = C'/ 


chiral state 


A^ = 


= -lKy = C, 4=4=0 


anisotropic state 


A, = 


= 4 = c. A., = 4 = 



Table 2: Examples of degenerate states of the renormalised mean field Hamilto- 
nian at half-filling, see also [20]. The general constant C = 0.479 is determined 
by the RMFT gap equations, ([Ml) and ([SU). *We note that the d-wave density 
matrix is not the d-density wave (DDW) order discussed in [5S] . 
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field ground states at half-filling and its possible relationship with competing or 
coexisting order are not yet fully understood. 

4.3 RMFT for the Hubbard model and application to 
HTSC 

The RMFT presented so far can be improved by considering all terms contribut- 
ing to the t-J Hamiltonian i.e. including the density-density term ^ as 
well as the correlated hopping terms pU]) . The inclusion of these terms is also 
necessary if one were to use RMFT for the Hubbard Hamiltonian, since the uni- 
tary transformation Ht-j = Pgc^^ He~^^ Pq (discussed in section between 
the Hubbard and t ~ J Hamiltonians lead to these terms. Physical quantities 
for the Hubbard model can be evaluated by considering expectation values in 
the re-transformed wave function e~^^ Pg\^o) |7j. Below we will use this ansatz 
to study the superconducting order parameter. In the following, we will include 
next nearest neighbour hopping {t') into the Hubbard Hamiltonian in order to 
allow for quantitative comparison with experimental data for the HTSC. 

4.3.1 Generalised gap equations for the strong coupling limit 

We obtain the RMFT gap equations for the Hubbard model in the strong cou- 
pling limit by considering the corresponding effective Hamiltonian, i. e. the full 
t-J Hamiltonian. The gap equations for this t-J Hamiltonian, which includes 
all terms from equation ([6]), can then be derived in the same way as described 
in the previous subsection. 

First we invoke the GA to obtain the rcnormaliscd Hamiltonian for ([6]). 
We note that all (nearest as well as further neighbours) hopping terms are 
renormalised by gt = (1 — n)/(l — n/2) and all superexchange terms by gs = 
1/(1 — n/2)^. Since the density-density term commutes with the projection 
operator Pq, it does not pick up any Gutzwillcr renormalisation factor. The 
new correlated hopping terms, equation (jlOp . are of the following form, 

(4+Ti,T4i'^^-L'^»+^2,T)PG*o ~ .93 (cj+^i,-f"-j,j,(l - '^»,t)cj+t2,T>*o > 

'\4+TuA,\(^^,i(^'>-+r2A)PG-'i'a ~ .93 (cj+^i44|Q4Q+^2,t)*o , (107) 

involve three sites (i, i 4- ti, and i + T2), and are renormalised by a factor 
33 = (1 — n)/(l — n/2)^. For a derivation of the GA for the correlated hopping 
terms we refer to the appendix of [S]. 

Next we decouple the resulting renormalised Hamiltonian by the same scheme 
discussed in the previous subsection, obtaining therefore the same gap equations, 
and as before. However, the dispersion relation, £'k = + A^, with 
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( i \ 

Ck = I 2,9tt + J -xi + J3—X2 I (cos kx + cos ky) 

( I' i \ 

— I 2gtt' + J' —xi + J3 -X2 j 2 cos k^ cos ky 

E ^ - ^2)] - M , (108) 

-ri#-r2 

Ak = J ^ [35,, + 1 - (3 + a::)53] (cos k^ - cos fcy) , (109) 

incorporates the effects of further neighbour hopping and correlated hopping 
terms. These expressions for and Ak are vahd for A = jA^^I = |Ay|, Aa; = 
—Ay, ^ = = i.e. for the d-wave pairing state, which is the most stable 
solution of the gap equations (|98|) and ((99)l . Note that we defined new hopping 
amplitudes for next nearest neighbours, ^' = ^x+y = £,x-y The last sum in pOSp 
is a sum over all pairs of non-identical neighbouring sites ti and T2, where t^-^ 
and are nearest and next nearest neighbour hopping terms. We defined, J = 
At^/U, J3 = At't/U, and J' = Af^/U and abbreviated, xi = 3^^ - 1 + 3(3- 2:)g3, 
X2 = 4(3 - x)g3, and x_d = (1 - x'^)g3 in (|108|) . 

As in section l4?2l the ground state jvE'o) of the renormalised t-J Hamiltonian 
results from above equations. By including a projection operator Pq into the 
wave function we obtain Pcj^t'o), which corresponds to a variational wave func- 
tion for the ground state of the t-J Hamiltonian in the fully projected Hilbert 
space. Invoking the canonical transformation e^*^ then provides an approxi- 
mate ground state e~*'^PG|^o) for the Hubbard model. 



4.3.2 Results from the generalised gap equations 

For a comparison with experiments, we follow [7] and consider a ratio t' /t = 
— 1/4 between next nearest and nearest neighbour hopping amplitudes, a value 
used widely in the modelling of the band structure of various HTSC [57] . Fur- 
thermore, we choose an on-site repulsion U = 12 1, i.e. we work in the strong 
coupling regime U ^ t,t', where the transformation from the Hubbard to the 
t-J model is valid approximately. The above choice of the model parameters 
reduces the number of free parameters to one energy scale, t w 300 — 500 meV, 
for the HTSC. 

In figure [29] we show the doping dependence of the superconducting gap, 
|Ak| at k = (tTjO), within RMFT, which resembles experimental observations 
quite well. However, the magnitude of the gap is overestimated by a factor of 
about 2 (see scaling factor a = 1/2 in figure within mean field theory. This 
overall mismatch is attributed to the fact that dynamical [155] and long-range 
correlations are neglected within RMFT, which is based on a local and static 
molecular-field approximation. 
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X 

Figure 29: Doping dependence of (solid line, right scale) the superconducting 
order parameter, $, and (dashed Hne, left scale) the superconducting gap, |Ak|, 
at k = (tt, 0) for t ~ 300 meV. The RMFT superconducting gap is scaled by a 
factor a = 1/2 for comparison with experimental data (red circles, Bi2122 [24]). 
From [7]. 

As mentioned in section I4.2.3| the superconducting gap is not identical to 
the true superconducting order parameter, $ = |(c||c|_,_^^ — c|^c|^_,_^|)| pOl [^ . 
Here, we determine the expectation value of $ within the re-transformed wave 
function, e~^^' PqI'^o} ■ Following section \TM we evaluate the canonical trans- 
formation e~''^ in order 0{t/U). Doing so, provides systematic t/f/-corrections 
to the result from section [4.2.31 where we used the wave function PcY^o) in 
calculating the expectation value of the superconducting order parameter. 

The calculations for the expectation value of a general observable O within 
the Hubbard model are summarised by, 

(O)e-.sp^*o = (e'^Oe-'^)p^*„ (110a) 
« (0+ z[5,6])p^*„ , (110b) 

where the last step corresponds to the evaluation of e~^^ in order 0{t/U); 
compare with section [2^ Note that ()110bp corresponds to an expectation value 
of the operator 0+ i[S, O] in the projected state PqI'^o) ■ We can therefore use a 
generalised Gutzwiller approximation by invoking the counting arguments given 
in section [01 

Setting, O ~ c|^|c|^_|_^| — c||c|_|_^|, the superconducting order parameter $ for 
the Hubbard model can be calculated using by above scheme. One finds, 

5,A + ^g3(6-x)A|. (Ill) 

In deriving (jllip we considered i' w within S, for simplicity, since \t'\ <C 
\t\. As shown in figure <!> vanishes as a; — ^ 0, and the i/J7-corrcctions do 
not qualitatively change the result of Zhang, et al. |2D] near half-filling. We 
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Figure 30: Doping dependence of Fermi velocity, vf- The RMFT results are 
compared with experimental data from Zhou, et al. [39], Borisenko, et al. |156j . 
and Kordyuk, et al. [157] , From [7]. 

emphasise that the above procedure can be used to calculate the expectation 
value of any observable and provides a systematic way to study the Hubbard 
model in the strong coupling limit. 

Next we consider the nature of the low lying excitations, i.e. the quasiparti- 
cles created at the nodal point, kp = kp^ = kp^y. The nodal dispersion around 
kp is characterised by the velocity, vp, which directly influences a number of 
experimentally accessible quantities. The Fermi velocity vp can be directly cal- 
culated by the gradient of along the direction, (0, 0) (tt, tt) within RMFT. 
The results obtained by using from equation (|108|) are presented in figure [30] 
(for t = 0.3,0.4,0.5 eV and ao = 4A). Edegger et al. [7] also found that the 
theoretically is well approximated by the formula. 



where J' and J' are set to zero for simplicity. As we will discuss later, the 
effective values of J' and J" within the GA become zero at half filling. So, 
ignoring their effect on the dispersion modifies the result only weakly. As seen 
in figure 1301 vp increases with x, but remains finite as a; ^ 0. In addition, 
we can infer from (|112p that the energy scale of the nodal velocity at a; = 
is determined by J, i.e. vp/{af)J) « ^/2sm{kp)^-£, « 1.5 (with ^ « 0.38 
and kp « The observed doping dependence stems from the effects of the 
Gutzwiller projection Pq. As x increases, holes gain kinetic energy by direct 
hopping, viz. gt increases with doping; but decreases, leading to the doping 
dependence of vp seen in figure [30] Note, that the RMFT results presented 
in figure [30] are absolute in value. No rescaling has been made for comparison 
with the experiments, contrary to the results for the gap |Ak=(7r.7r)| presented 



vp/ao ~ \/2smkp 2gt{t + 2t' coskp) + xi—S, 



(112) 
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Figure 31: Energy dispersion, oj = — i?k = -\/^k + Ak, of the Gutzwiller- 
Bogoliubov quasiparticle for different doping x. From [7]. 

in figure [29l 

The above resuhs are in agreement with the VMC resuhs of Paramekanti, et 
al., who extracted vp from the discontinuity of the first moment of the spectral 
function in the repulsive U Hubbard model ([2T | 1117] . also section [5.3.2|) and 
of Yunoki, et al., who obtain vp from the quasiparticle dispersion in the t- 
J model ( 11461 . also section \T72\ . These results also yield a good fit to the 
ARPES data [3i [TSl [TS71 [UHl [HSl , as illustrated in figure [301 We note that 
the doping dependence oi vp in the severely underdoped regime remains to 
be settled experimentally. While some groups report a nearly constant Fermi 
velocity (see data for LSCO in figure [3D]), others observe a slight increase with 
doping (see data for YBCO and Bi2212 in figure [30]) . We further emphasise 
that the energy scales t and J might be extracted from the ARPES data in 
Vp, whenever data with high accuracy becomes available. By using ^ « 0.38 
and setting kp and the ratio t' /t to the experimentally observed values, t and J 
can be fitted by (|112p . In addition, the RMFT calculations find that the nodal 
properties remain essentially unchanged when A is set to 0, i.e. the doping 
dependence of vp results from the vicinity of the RVB state to a Mott insulator, 
rather than the occurrence of superconductivity itself. 

In figure [3T1 we present the energy dispersion, lj = — i?k: of the Gutzwiller- 
Bogoliubov QP along the directions, (0,0) (tt, 0), (7r,0) — *■ (7r,7r), and 
(tt, tt) (0, 0) for different doping levels x. The dispersion is flattened when 
approaching half-filling and the gap around (tt, 0) becomes large. We empha- 
sise that these RMFT calculations adequately describe only the low energy 
sector of the HTSC, and do not seek to explain the 'kink' at higher energies 
[39l [T56l[T57l[T58l[T59] . 

Equations (|108p and (|109p can also be used to study the effects of t' onto 
the magnitude of the superconducting order parameter Figure [32] shows 
the value of $ at optimal doping as a function of the ratio t' /t in the bare non- 
interacting dispersion. We observe a maximum at about t' /t = —0.35 depending 
slightly on whether we seiF^ t" = or t" = —t'/2. This observation is in 

^■'We note that the inclusion of t" = —t'/2 into the bare dispersion is sometimes used to 
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Figure 32: The magnitude of the optimal superconducting order parameter (opt. 
$) as a function of the ratio t' /t as obtained in the RMFT calculations for the 
Hubbard model. Opt. $ provides the value of <i> at optimal doping, i.e. the 
maximal value of $ for this particular set of model parameters (t'/i, U ~ 12t). 
The dashed line shows calculations, where third nearest neighbour hopping has 
been taken into account into dispersion by setting t" ~ — i'/2. 



agreement with band-structure calculations [160] where it is found empirically 
that compounds with ratio t' /t « —0.1 in the dispersion (determined from 
band-structure calculations) have smaller Tc (corresponds to smaller <&) than 
compounds with a ratio t' /t « —0.35. The RMFT calculations in figure [35] also 
match VMC results of Shih et al. |161| , which we will discuss in more detail in 
section 15.2.21 

To summarise, the calculational scheme we described above, presents a sys- 
tematic way to study the Hubbard model in the strong coupling limit. It is 
based on the validity of the t-J model as an effective Hamiltonian (in the large 
U limit) and on determining expectation values within the re-transformed trial 
wave function, e~''^PG|^o); a scheme that can be extended to excited states as 
will be described in section [6] (see also [7]). 

4.4 Possible extensions and further applications 
4.4.1 Incorporation of antiferromagnetism 

The incorporation of antiferromagnetism is an example of a possible extension of 
the RMFT. In order to describe an antiferromagnet with finite sublattice mag- 
netisation m, we have to allow for an additional degree of freedom in the wave 
function. When deriving the corresponding gap equation we must keep in mind 
that the antiferromagnetic correlation affects the GA as discussed in section [XT] 
However, Himeda and Ogata [102] showed by VMC calculations that even the 
formulas from section 13.11 do not adequately describe all aspects in a magnetic 

get a better fit with band-structure calculations and ARPES data. 
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Figure 33: RMFT calculations including antiferromagnetism from Ogata and 
Himeda |107j . The self-consistent parameters A (pairing amplitude), ^ (hopping 
amplitude), and m (staggered magnetisation) arc shown as a function of the 
doping rate a; = 1 — n for J /t = 0.3 and t' = J' = 0. The dashed lines represent 
the results when the antiferromagnetic order is suppressed, i.e, m is fixed to 
zero. From |107j . 



ordered state. They determined effective Gutzwiller renormalisation factors by 
comparing the numerically obtained expectation values in the projected state 
with the respective mean field values before projection. It was found that the 
z-component of the Gutzwiller renormalisation factor is enhanced compared 
with those of the xy-component . 

Ogata and Himeda |107j argued that the discrepancies stem from spatial 
correlations neglected by the GA. They derived extended Gutzwiller renormal- 
isation factors by considering a cluster around the sites i and j to incorporate 
further inter-site correlations. Applying these renormalisation factors and solv- 
ing the gap equations including antiferromagnetism yields the results of figure 
[33l We see that for doping, 5 < 0.1, long range antiferromagnetic order coexists 
with superconductivity. For higher doping the magnetisation m vanishes and 
solely the superconducting order remains. This result is obtained neglecting the 
next nearest neighbour hopping {t' = 0) and agrees with previous VMC results 
[1021 11031 11041 1105] . We note that the extended Gutzwiller renormalisation fac- 
tors of Ogata and Himeda are essential for reproducing the VMC calculations. 
However, figure [33] does not quantitatively agree with the experimentally phase 
diagram of hole-doped Cuprates, where antiferromagnetism disappears at about 
3 — 5% doping. A better match may be obtained by considering the effects of t' 
as done in VMC [106] and quantum cluster studies [1621 1163] . 
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4.4.2 Applications to inhomogeneous systems 

The RMFT described has also been used to study inhomogeneous phases such as 
stripes and checkerboard charge order [SSI IMl HOOI I101| , vortex cores |llll 1112] , 
magnetic and non-magnetic impurities [1081 11091 IllOj . These investigations 
throw light on the interplay between antiferromagnetic correlations, d-wave su- 
perconductivity, and charge order and can be compared with STM data. 

However, such studies require an unrestricted Hartree-Fock treatment of 
the rcnormaliscd t-J Hamiltonian (j93p . i.e. local expectation values such as, 
Ay = (c^|C^j2*o: s,TLd, ^ija = {cl^Cjc)\si„ , must be considered independently for 
each bond [164| . Furthermore, the local charge densities riicr generally differ from 
site to site, and thus the Gutzwiller renormalisation factors of the renormalised 
Hamiltonian depend on the site indices i and j {g\^ , 55 ) • Special attention 
must be paid when deriving these Gutzwiller renormalisation factors, because 
the local charge densities can differ between the projected and unprojected state 
(see discussion in section 13. ip . For inhomogeneous systems the RMFT gap 
equations generalise to the so-called Bogoliubov-dc Gcnncs equations, which 
must then be solved self-consistently. 

The investigation of charge modulations within above framework [98 t [99 lfTom 
llOlj provides an understanding of the 4x4 checkerboard patterns seen in the 
STM data of the HTSC. These studies neglect long range antiferromagnetism 
and assume S^ij = = S^iji and rti^ = riij. This is a reasonable assumption 
since the authors concentrated on doping levels where antiferromagnetism is not 
observed experimentally. The renormalised mean field Hamiltonian can then be 
written as [lOOj . 

Hmf = - ^ X! 9t (.c\,aCj,<j + h.c.) - ^ n,^a 

{■ij)a la 

- J^E5l(e,4.c,,. + h.c.-ii,,p) 

- JjE.g|(A,4,c],_,+h.c.-|A,,f). (113) 

However, we must abandon above constraints, = ~ ^iji and = 
riii, for investigations around vortex cores or impurities, where antiferromag- 
netic correlations are essential. Doing so and solving the Bogoliubov-de Gennes 
equations under an uniformly applied magnetic field shows that significant an- 
tiferromagnetic correlations develop inside vortex cores jll2| in agreement with 
experimental observations [HSl [HSl [Ml MB \lM [HO]. Tsuchiura, et al. 
[108. 1091 1110] also studied the effects of magnetic and non-magnetic impurities 
onto the local density of states in HTSC within above approach. The results 
obtained resemble the STM data [1S1[171[1H| quite wel0 The self-consistent de- 
termination of order parameters within the renormalised Bogoliubov-dc Gennes 

^^Nevertheless VMC calculations 11711 for a non-magnetic impurity report some minor dis- 
crepancies to the corresponding RMFT study I110| . 
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theory was also applied to study surface effects in 2D superconducting states 

[niiizs]. 

To analyse the above problems within an unrestricted Hartree-Fock theory, 
most authors consider a large (but finite) unit cell, which exhibits a certain 
charge ordering pattern or which has a vortex core or an impurity site in the 
middle. The corresponding renormalised Bogoliubov-de Gennes equations can 
then be solved by assuming a lattice of unit cells {e.g. iVc = 20 x 20) and 
making use of Fourier transformations. While most studies use the Gutzwiller 
factors derived in section [3] some recent works [1121 1110] use the extensions 
proposed by Ogata and Himeda [107j (see section 14.4. ip . Finally, we note that 
all these studies concentrate on the ground state properties (T = 0). It would 
be very interesting to consider finite temperature effects within a renormalised 
unrestricted Hartree-Fock theory, and to our knowledge, such studies have not 
yet been carried out. 

4.4.3 Gossamer superconductivity 

Another class of renormalised mean field theories considers a modified version 
of the Hubbard model, which includes a superexchange interaction J like in the 
t-J Hamiltonian. This t-J-U model was proposed by Zhang |174j to study the 
so-called gossamer superconductivity [175] . Here, the form of the GA, which 
includes finite double occupancy, must be used for the renormalised Hamiltonian 
|174j . The RMFT gap equations are obtained in a straightforward manner, 
and the number of double occupancies is determined by optimising the ground 
state energy. Within this approach, at half-filling, there is a first order phase 
transition from a Mott insulating phase at large Coulomb repulsion U to a 
gossamer superconducting phase at small U. Away from half-filling the Mott 
insulator evolves into an RVB state, which is adiabatically connected to the 
gossamer superconductor |176| . Some authors follow this approach to study 
HTSC [1761 [ml [ml [Hi [164] while others have used it in the phenomenology 
of organic superconductors [HH [23] ■ 

4.4.4 Time-dependent Gutzwiller approximation 

The studies discussed so far mainly focused on the superconducting state. Sei- 
bold and Lorenzana |180| considered the Hubbard model without superconduct- 
ing pairing and developed a time-dependent Gutzwiller approximation anal- 
ogous to the time-dependent Hartree-Fock theory. This formalism incorpo- 
rates ground state correlations beyond the GA within the random phase ap- 
proximation and allows for a computation of the dynamical density-density re- 
sponse function. The scheme successfully describes several interesting features 
of HTSC, such as the dynamics of stripes [181| or the dispersion of magnetic 
excitations [1821 11831 1184j and was recently applied to investigate checkerboard 
inhomogeneities |185j . It would be very useful if this scheme could be adapted 
to study the dynamics of a homogeneous superconducting phase. 
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5 Variational Monte Carlo calculations for HTSC 
- an overview 



The VMC technique aUows for an accurate evaluation of expectation values in 
Gutzwiller projected wave functions. In this section we present technical details 
of the VMC method and review the variational search for the optimal ground 
state energy in the Hubbard and t-J model. In this context, we also discuss 
the coexistence of superconductivity with antiferromagnetism and flux states 
as well as the improvement of the trial wave function by further variational pa- 
rameters and Jastrow factors. Further we consider doping dependent features of 
projected wave functions and compare them to experimental findings in HTSC. 
Finally, we discuss a recent numerical study dealing with the tendency towards 
a spontaneous breaking of the Fermi surface symmetry. 

5.1 Details of the VMC method 

The VMC method was first applied to the study the projected Fermi sea |116i 
1186] . which has a fixed particle number. However, superconducting BCS wave 
functions |^'o) are generally defined in a grand canonical ensemble, where the 
wave function shows particle number fluctuations. These particle number fluc- 
tuations are also present in the projected BCS wave function, \^) = Pq |^o)- 
Within the VMC scheme we have now two possibilities, we can either directly 
work with \ ^) (grand canonical ensemble) [97], or we can project out the particle 
number fluctuations by a projector Pn that fixes the particle number and work 
with l^*) = Pn\^) (canonical ensemble) [17]. In this review, we only present 
the method of Gros [T7| , I147j . used in most recent VMC calculations, since it 
avoids complications due to the fluctuating particle number. The possible dis- 
crepancies between the grand canonical and the canonical VMC scheme have 
been discussed in detail in section [3?2l 

5.1.1 Real space representation of the trial wave function 

Before performing a VMC calculation we have to rewrite the wave function in 
an appropriate way. By inserting a complete set of states { |q;) } for the subspace 
that excludes double occupancy (and with a fixed particle number N), we can 
remove the projection operator P^ and Pq from the wave function j^P), 

1^-) = FatPgI^o) = 5I("l*o) \a) ■ (114) 

a 

The most suitable choice of \a) is given by a straightforward real space repre- 
sentation in terms of fermion creation operators. 

The state (|115p is specified by two disjoint sets {Ri . . . Rjv^ } and {Ri . . . Rat^ }, 
which determine the positions of the up- and down-spin electrons on a finite 
lattice. 
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Next we have to calculate the overlap (a|^'o). To determine this quantity by 
a Monte Carlo calculation, we write the BCS wave function P^vl^'o) as [15], 



PatI^o) = PNl[[uk + Vkcl^cl^^j\0) (116a) 

k 

« PNl[(l + akcl^cl^^\0) (116b) 



N/2 

cx ( >:akcT^cLk| I |0) (116c) 

\ N/2 

J2 a(R,,i-%T)4,,„T4,,,a |0) . (116d) 
In (|116b[) we defined the quantity Ok = Wk/uk, which can be written as 



A' 



using the mean field result from ([W]) . Due to the projection operator Pjv we 
can then represent the wave function by a product of N/2 pairs, where we use 
No- ^ = N[ = N/2, valid for a BCS wave function. In (|116dp we assumed 
ak — fl-k, applied a Fourier transformation, and defined 

a(r) EE ^ flk cos(k • r) . (118) 

k 

Finally, we arrive at the real space representation of PnI'^o) as in (|116d[) . 

Since all configurations a in (|114p have the same to particle number N, 
(a|*o) = (a|Piv|^'o)- Making use of (|116dp one finds that the overlap, (a|*o), 
is given [171 1147j by the determinant of the matrix which has the form, 

f a(Ri,x-Ri.T) a(Ria-K2,T) ■■• a(Ria - R^v^.t) \ 
a(R2,x-Ri,T) a(R2,i ^ R2,t) a(R2a ^ Ra^,,,!) 

\ a(Rw<,a ~ Ri,t) a(Rw„,j, - R2,t) ■•• a(Rw,,a ~ Rjv^,t) / 

To see this we must expand (|116dp and gather all terms contributing to the 
configuration a, which has down-electrons on {Ri |, R24 . . . Rat^j} and up- 
electrons on {Ri_f,R2.| . . .Rat^^^i}. The number and functional form of these 
terms are obviously the same as those for \Aa\. Next we must order up- and 
down-electrons in the same way for all terms. By doing so we pick up relative 
signs, which are exactly reproduced by the determinant of Aa. 

We note that the above real space representation can be extended [86] to 
wave functions, which allow for a staggered magnetisation and an unequal num- 
ber of up- and down-electrons, N^ 7^ ^i- Then, the a(r) in (jll6d[) and pi8p 
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becomes spin and site dependent, i.e. a(r) — > a(Ri.£ri , Rj.o-j , Ci, CTj). The val- 
ues of a(Ri_CTi J R-j.cTj ) CTi, CTj) depend on the particular choice of the mean field 
wave function and can be evaluated numerically. The overlap (aj^'o) is then 
determined by [5S] . 

(alfo) - P/(Q) , (119) 
where P/(Q) is the PfaffiarF^ of the matrix 

Q = a(Ri^cri,Rj,o-j,o-i,o-j) - a(Rj>j,Ri,cr,,(Tj,(Ti) . (120) 

The positions of the electrons, Ri,o-i and Rj,<Tj determine the real space configu- 
ration a. For a simple BCS wave function with a(Ri_|, Rj.t, Ti T) = i 
, i) = and a(Ri_| — R-j.i) = ^(R'j.Ti F^i.i: Ti l)i the overlap (a|\I'o) in equation 
(|119[) reduces to the previously discussed determinant \Aa\. 



5.1.2 Implementation of the Monte Carlo simulation 

Using (|114p . we may write the expectation value of an operator O in j'l') as, 

(6). = (^fPcP.dP.Po\^o) 



y^^^mpP^^^ (121b) 



[y^ (a|d|/3)(/3|^o) \ \{^^\a)? 
^ \^ (vfola) J (*o|PjvPg|*o) ^ ^ 

= ^/(a)p(a), (121d) 



with. 



= , 'f^°'"^f , . (122b) 
Here, a and /3 are real space configurations (|115[) . Since, 

p(a)>0, ^p(a) = l, (123) 

a 

are the features of a probability distribution, we can evaluate (O)* by a ran- 
dom walk through the configuration space with weight p(a). Therefore, we can 



^^The Pfaffian is the analogue of a determinant which is defined only for antisymmetric 
matrices. For an antisymmetric matrix A, the square of the Pfaffian is equivalent to its 
determinant, viz. Pf{A)^ = 1^1- 
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Figure 34: (left) Real space picture of the + 1 lattice for L ~ 5, with periodic 
boundary conditions, (5, 1) and (—1, 5), apphed along the opposite edges of the 
tilted square indicated by dashed lines, (right) The k-space Brillouin zone of 
the 'tilted lattice' for L = 5. From pTZj . 

analyse (|12ip by a standard Metropolis Monte Carlo calculation. We note that 
the norm (4'o|-P/v-Pg|^o) in (|122bp is not of relevance within the Monte Carlo 
calculation, since only relative probabilities p(a) enter the transition probability. 

Next, we make comments on the updating procedure and the calculation 
of the determinant \Aa\. Most VMC calculations generate a new configuration 
a' by randomly interchanging two electrons with opposite spin or moving an 
electron to an empty site. The random walks thus constructed are ergodic. 
In general, to optimise the numerical performance, the rules for generating the 
random walk through the configuration space should be chosen so as to maximise 
the acceptance rate, T{a a'). 

The calculation of \Aa\ is numerically expensive and is required at each 
Monte Carlo step for the computation of p{a). Therefore it is advantageous to 
determine the ratio between new and old determinant (new and old 

configuration a' and a) instead of directly evaluating \Aai\ for every configura- 
tion. According to Ceperley, et al. [1 87] . this ratio can be efficiently computed 
within 0{N^) computation steps, while a direct evaluation of \Aar\ requires 
0{N^) steps. The trick is to store not only the matrix Aa, but also its inverse 
A~^. For the commonly used updating procedures mentioned above, a' differs 
from a only by the interchange of two electrons with opposite spins or the inter- 
change of an electron and an empty site. Thus, the matrices Aa> and A^ differ 
only by one row and one column, and |j4q/|/|Aq| — \AaiA~^\ , which enters the 
transition rate T{a a'), can be easily computed. 

A general advantage of the Monte Carlo method is the possibility to estimate 
the numerical accuracy systematically with the error being proportional to the 
inverse square root of the number of Monte Carlo steps Nr. Present computer 
capacities allow us to consider sufficiently large clusters, where finite size effects 
play a minor role. However, Ok = vu/uk, defined in equation (|117p becomes 
singular whenever Ak = and ^k < 0. In particular, this is problematic for a 
c?-wave order parameter, for which Ak = for all A:-points along the Brillouin 
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zone diagonals, i.e. \kx\ = \ky\. It is thus convenient to avoid these k-points 
by an appropriate choice of boundary conditions. There are three different 
approaches discussed in hterature. The first possibihty is a tihed lattice with 
periodic boundary conditions (PBCs). Such a lattice has + 1 sites with 
odd L, preserves the fourfold rotational symmetry of the lattice, and does not 
introduce any twist in the boundary conditions. An example of these (widely 
used) boundary conditions (see, e.g. [r7llll7[[2Tlll47j ) is illustrated in figure [34l 
Another choice is an L x L lattice with even L and periodic and antiperiodic 
boundary conditions in the x- and the y-direction, respectively. Finally it is 
possible to use a rectangular x Ly lattice with PBCs and mutually coprime 
dimensions and Ly, i.e. the greatest common divisor of and Ly being 1. 

5.2 Improvements of the trial wave function 

The early VMC calculations for projected ECS states of Gros [T7| and Yokoyama 
and Shiba |97| were carried out to check whether a Gutzwillcr projected super- 
conducting wave function is indeed a variationally good trial state for the t-J 
model. To limit the number of variational parameters, these authors used a 
dimensionless dispersion, = —2 (cosfc^^ + cosfcj^) — /.t, and various supercon- 
ducting gap functions Ak to calculate ak = Ak/(^k + \/^k~+~^)- -'^'^ orig- 
inal work, Gros |17| compared variational energies of s-wavc, Ak = A, rf-wave, 
Ak = A (cos fca; — cos fey), and extended s-wave, Ak = A(cosA:2; -I- cosfcy) — fi, 
functions. By optimising solelj{l3 the variational parameter A, he found that a 
d-wave gap can substantially lower the energy compared to projected Fermi sea 
(Ak = 0) at half-filling as well as at finite doping. The result is consistent with 
other early works such as the VMC calculations of Yokoyama and Shiba [97] or 
the mean field theories of Zhang, et al. [20] and Kotliar and Liu [18]. 

More detailed studies showed that the optimal superconducting state re- 
mains a pure d-wave even when mixed states of s- and d-wave pairing are con- 
sidered [188] . The optimal variational parameter A decreases when going away 
from half-filling and vanishes at about 30% doping. The exact dimension of the 
superconducting region in the phase diagram depends on the choice of J/t as 
well as on the inclusion of the correlated hopping term [188j . 

5.2.1 Antiferromagnetism and flux states 

Further extensions [102111031 11041 [TU51ll06j considered the incorporation of anti- 
ferromagnetism for a more accurate description of the t-J model near half-filling. 
These studies show a coexistence between superconductivity and antiferromag- 
netic long-range order (AFLRO) for doping x < 0.1. At half-filling, the optimal 
so-constructed wave function has a staggered magnetisation of 0.75 and a varia- 
tional energy of —0.664 J per site; impressively close to the best numerical esti- 
mate of —0.669 J per site by Green's function Monte Carlo techniques [18911190] . 

^^Due to the fixed particle number, the chemical potential fi becomes an additional free 
parameter. However, this parameter was fixed in |17| by setting the chemical potential fi to 
those of the unprojcctcd wave function. 
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Figure 35: Optimal energies per site (in units of J) for five different Gutzwiller 
projected wave functions with a linear part subtracted /iscp) as a function of 
doping X. Wave functions: superconducting without antifcrromagnetism (SC, 
solid circles), superconducting with antifcrromagnetism (SC+AF. empty cir- 
cles), staggered-flux without antiferromagnetism (SF, soHd squares), staggered- 
flux with antiferromagnetism (SF+AF, empty squares), and zero- flux (projected 
Fermi sea, empty diamonds). The arrow in the panel shows the best variational 
estimate for the half-filled system {E = -0.669 J per site) [HSl [190]. Only 
nearest neighbour hopping is considered and J/t ~ 0.3. From |105j . 
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Figure 36: Momentum distribution n(k) in the first Brillouin zone for doping 
X = (a) 0.19, (b) 0.31 (c) 0.42 (d) 0.49 for 12 x 12 lattice with J jt = 0.3, 
t' = — 0.3<, and i" = 0. (e) optimal parameters t'^ (squares) and A (circles). 
From [T6T] . 

A comparison of the variational energies of the different wave functions is given 
in figure 1351 The figure also reveals an upward convexity of the ground state 
energy (SC+AF state) as a function of doping. This indicates a phase sep- 
aration at ccscp = 0.13 (see figure [SH). However, Ivanov [105j demonstrated 
that a sufficiently strong nearest neighbour Coulomb repulsion can suppress the 
formation of separated phases. Further VMC calculations showed that the co- 
existence of superconductivity and AFLRO is nearly absence if next and second 
nearest neighbour hopping are included [106j . For these more realistic model 
parameters, the AFLRO disappears at about 6% doping in better agreement 
with experimental observations [106j . 

Apart from the superconducting states, the projected staggered-fiux state 
has also been studied as a competitive variational state; however, its energy lies 
above those of the c?-wave for all dopings (figure [35|). As discussed in section 
m the flux state becomes identical to the superconducting state at half-filling 
explaining the collapse of the energies in figure [35] (see also [142j ). This be- 
haviour is due to S't/(2)-symmetry, which is also responsible for the occurrence 
of staggered- vorticity correlations of current in the d-wave state at small dopings 

m- 

5.2.2 Increasing the number of variational parameters 

In recent VMC calculations, the chemical potential fi as well as the next nearest 
neighbour hopping t'^^^j. are chosen as additional variational parameters, which 
are optimised numerically. While the chemical potential has minor influence 
on the optimal state |188j . a variational t'^^^^. can signiflcantly effect the shape 
of the Fermi surface. Himeda and Ogata [192] reported that for a bare disper- 
sion t' ~ and a doping level of x = 0.15 the lowest energy is provided by a 
variational t[,^^ = ~0.lt, causing a spontaneous deformation of the Fermi sur- 
face. More detailed VMC studies [161] include next nearest (t') and next next 
nearest (t") neighbour hopping in the bare dispersion and also use variational 
parameters t[j^^ and t'^i^j.. The obtained momentum distribution n(k) (related 
to the Fermi surface, see section [^3]) together with the optimal variational A 
and t(,a^j from these calculations are illustrated in figure [36l This work of Shih, 
et al. |161j also revealed that the more negative the bare ratio t'/t, the higher 
the superconducting pairing in the optimal variational state of the t-J model. 
This is in agreement with band-structure calculations |160) that suggest that 
the ratio t' /t is essential to raise Tc- Similar trends can be inferred from the 
RMFT calculations for the Hubbard model discussed in section [4.3.21 

VMC studies of inhomogeneous phases |193| find that around x = 1/8, 
stripe states with fluctuating c?-wave superconductivity can lower the variational 
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Figure 37: Fermi surface of the isotropic t-J model with J = 0.3t and t' = —0.3t 
at a; = 0.07 (a) RMFT results for the Fermi surface of the normal state with 
Ak = (quasi ID state, solid line) and the optimal d-wave state (isotropic, 
dashed line), (b) Best quasi ID state on a (15^ + l)-sitcs lattice by VMC; filled 
circles indicate the Fermi surface. From [8|. 



ground state energy in the two-dimensional t-J model. More recent studies 
report that at x = 1/8, a bond-order modulated staggered flux state can also 
overcome the RVB superconductor for sufficiently large short range Coulomb 
repulsion [194j . However, the energy gains within these studies arc often quite 
small and sensitively depend on model parameters. Nevertheless, these VMC 
calculations show that the slightly doped t-J model exhibits tendencies towards 
various inhomogencities, which could be relevant for explaining some of the 
experimental observations in the underdoped HTSC. 

The energy of the projected d-wave state can be improved further by the 
incorporation of Jastrow factors (see section 12. 6p . Sorella, et al. |145| showed 
numerically that such wave functions lower the variational energy and still ex- 
hibit long range superconducting order. Nevertheless there is still debate (see 
|1951 1196] I197j ) whether the superconductivity within the VMC scheme results 
only from a biased choice of the wave function or is indeed a ground state prop- 
erty of the t-J model. In our opinion, this debate docs not pose an obstacle 
towards our understanding of the HTSC, since we are primarily interested in 
physical properties of projected wave functions rather than in proving them to 
be exact ground states of a particular Hamiltonian. In this, we follow the point 
of view espoused by Anderson, et al. [198| . 'The philosophy of this method is 
analogous to that used by BCS for superconductivity, and by Laughlin for the 
fractional quantum Hall effect: simply guess a wave function. Is there any better 
way to solve a non-perturbative many-body problem?' 

5.2.3 Investigation of the Pomeranchuk instability 

The possibility that strong correlations may break the symmetry of the under- 
lying Fermi surface was studied recently in [S] . As illustrated in figure EUa) , 
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Figure 38: (a) VMC results for condensation energies per site Ccond of the 
quasi ID state ((5^° 7^ 0, Ak = 0) and the d-wave state (Ak 7^ 0, (5^° = 0) 
with t' = -0.3 i; see and (fHS)) for the definition of these states. The 

optimal variational A^fr of the d-wave is shown in (b), the optimal variational 
asymmetry 5}^^ of the quasi ID state is given in (c). The errors in (b) and (c) are 
AAifi = 0.05 and AfJi^ ^ o.05, respectively. System sizes: L = 11^ + 1 = 122, 
L = 132 + 1 = 170, L = 152 + 1 = 226, and L = 17^ + 1 = 290. From g]. 



this instability results in a deformation of the Fermi surface, which becomes 
quasi one dimensional, although the underlying two dimensional (2D) lattice 
is still isotropic. Motivated by the Fermi surface depicted in figure ISTT a). the 
state resulting from the tetragonal symmetry breaking can be called a 'quasi ID 
state'. This phenomenon is also called a Pomeranchuk instability of the Fermi 
surface. 

To investigate instabilities towards quasi ID states, we have to extend the 
variational space by an additional order parameter, which allows for a finite 
asymmetry in the wave function. A possible choice was proposed by Edegger et 
al. [Hj, who determined Ok = Ak/(^k + V^Ck+^k) by using 

Ck = - 2 [(1 + $1^^) cos fc^ + (1 - 5l^,) cos ky] 

- At'^^^ cos cos ky - ^var (124) 

and 

Ak = A^f|.(cos k^ - cos ky) + ^iti^os k^ + cos ky) . (125) 

The five variational parameters, viz. the asymmetry (5^^ between the x— and the 
y— direction, the variational next nearest neighbour hopping term t'^^^, a varia- 
tional chemical potential /Xvar, and variational parameters for d— and s— wave 
pairing, Avar and Avar, can be optimised by determining the energy expec- 
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Figure 39: RMFT and VMC results for the d + s-wave ground state of the 
anisotropic t-J model with J = 0.3t and Sq^ = [t.^ — ty)/{tx + ty) = 0.05. (a) 
Effective asymmetry = {i^ — iy)/{ix + iy) from RMFT as a function of hole 
doping X for (dashed) t' = and (solid) t' = — 0.3i. VMC results for t' = are 
given by squares and circles for L — 122 and L = 170, respectively, (b) RMFT 
Fermi surface (solid lines) of the d + s-wave ground state and the tight binding 
dispersion (dashed) at a; = 0.08 with t' = -OM and ^J^d = 0-025. From [8]. 

tation values for different choices of variational parameters within a standard 
VMC technique. 

VMC calculations for the isotropic t-J model [8] show in agreement with 
previous studies [199(l200tl201j . that the optimal variational state remains a pure 
d-wave without any anisotropy (see figure [38]). However, when restricting solely 

to non-superconducting solutions, i.e. setting Ak = by using A^fr = A^lr = 0, 
a projected anisotropic Fermi sea provides a better energy than the isotropic 
one. In figure [551 this effect is shown for a VMC calculation [Hj in an isotropic 
t-J model. The figure also illustrates that the optimal d-wave state has much 
better energy than the quasi ID state, which is the best state in the variational 
subspace Ak = 0. 

The situation can be quite different when the underlying lattice structure 
is anisotropic. In this case, the tendency towards a quasi ID state is present 
even in the superconducting state. RMFT, VMC [8], and SBMFT pM [200] 
calculations predict an optimal state in which the bare anisotropy i5q^ of the 
lattice can be significantly enhanced due to the electron correlations. As seen 
in figure [Sgf a) . the bare asymmetry of Sq^ = 0.05 increases within the RMFT 
calculations up to about 5^^^. ~ 0.2 in the underdoped regime. These results are 
confirmed to some extent by VMC calculations [figure [5W a). circles and squares], 
that show an increase of the asymmetry up to about (5^^ w 0.1. Furthermore, 
the enhancement of anisotropy may even lead to a change in the topology of 
the underlying Fermi surface as can be inferred from figure [SOl b) . 

The PI is one out of several possible instabilities in the t-J model that arise 
from the effects of superexchange and that can be revealed by VMC and RMFT 
techniques. Since J oc /U , a small asymmetry in the bare hopping integral 
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Figure 40: Doping dependence of the superconducting order parameter $ from 
VMC calculations of Paramekanti, et al. [21]. Model parameters: U = 12t, 



t becomes twice as large in the superexchange energy. Hence, it is natural that 
the effects discussed in this subsection are largest in the underdopcd regime, 
where the dispersion is mainly determined by J. 

5.3 Ground state properties - VMC results 

Within this subsection, we discuss ground state properties of the HTSC by con- 
sidering observables in a Gutzwiller projected superconducting state. We follow, 
in part, Paramekanti, et al. |117i [21], who studied the Hubbard model in the 
strong coupling limit using the re-transformed trial wave function, e^^^ PqI'^o) 
(see section Hj). By evaluating the canonical transformation e"''^ to 0{t/U), 
this ansatz can be treated within the VMC scheme. The i/J7-corrections due 
to e~*'^ provide a more accurate description of the HTSC, however, the quali- 
tative nature of the results is not changed compared to the t-J model. In the 
following, we ignore the possibility of the superconducting state coexisting with 
a flux state, antifcrromagnetism, or a charge ordered state. 

5.3.1 Superconducting gap and order parameter 

In the previous subsection, we saw that the variational parameter A, which 
is proportional to the superconducting gap Ak, increases in the limit of half- 
filling. The doping dependence of A is illustrated in figure and resembles the 
RMFT result (figure [22]) as well as the experimental observed gap at k = (tt, 0). 
However, we cannot deduce the relevant energy scale of the gap from A, since 
it is a dimensionless parameter within the VMC calculations. For detailed 
statements about the gap we have to consider the energy of excited states as we 
will do in sections [6] and [7] 



t' = -t/4. 
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Figure 41: (a) The momentum distribution function n(k) plotted along the 
diagonal k = (fc, k) showing the jump at kp which implies a gapless nodal 
quasiparticle with spectral weight Z. (b) Nodal quasiparticle weight Z{x) as 
a function of doping x compared with the simple SBMFT result Z^^{x) ~ x. 
Model parameters: U = 12t, t' = -t/A. From pT] , 



When considering the variational parameter A, we must realise as discussed 
in section 01 that it docs not correspond to the true superconducting order 
$ = |(c-|C-_|_^j^ — c^cj_,_^|)|. The relevant physical quantity here is the off- 
diagonal long range order (ODLRO) [T71 [Wl ^ defined by, 

i^a,^(r-r') = (c|(r)c|(r + d)c|(r)cT(r + /3)) , 

where a,(3 = x,y. In the limit of large |r — r'|, Fa^p is related to $^ via 
Fa^p ±<i>^ with + (— ) sign obtained for d || (_L) to b, indicating c?-wave 
superconductivity [21]. The doping dependence of the superconducting order 
parameter (f> is depicted in figure HDl (VMC calculations of Paramekanti, et al. 
PT]V It is not identical to A as first noted by Gros [13 [UJ. The VMC 
calculation match the RMFT result (figure where $ vanishes linearly as 
a; — > 0. The vanishing order parameter $ indicates a Mott insulating phase 
at a; = 0, where superconductivity is destroyed by the suppression of particle 
number fluctuations. At finite doping x a superconducting state is realised in 
the range < x < 0.35. 

5.3.2 Derivation of spectral features from ground state properties 

Next we follow [21] and analyse the onc-particlc spectral function j4(k, lo) by 
calculating the moments, 

Mi{k) = dLuJA{k,Lu), (126) 



81 



in the projected d-wave ground state at T = 0. This ansatz allows to obtain 
information about A(k, w) from ground state expectation values without the 
need for explicit representations of the excited states. We first concentrate on 
the zeroth moment Mo(k) = n(k), which is equivalent to the moment distribu- 
tion n(k). Figure BTT a) shows that n(k) has a jump along (0,0) to (tTjTt). This 
implies the existence of gaplcss quasiparticles and allows us to write the low 
energy part of the spectral function along the diagonal as, 

A{k,Lo) = Z5{Lu-^k)+A'''\ (127) 

where, — vp{k — kp), is the quasiparticlc dispersion and A™'^, a smooth 
incoherent part. The location of the discontinuity determines the Fermi point kp 
and its magnitude, the quasiparticlc weight, Z. Figure HU b) reveals a significant 
doping dependence and shows that Z vanishes when approaching the Mott- 
Hubbard insulator x = 0. This behaviour is in agreement with more direct 
calculations, which explicitly include quasiparticlc states (sections [6] and [7]) as 
well as experiments. 

To determine the nodal Fermi velocity vp, we have to evaluate the first 
moment Afi(k) = {cl^^[H , ci^a-]) along the nodal direction. Due to the singular 
behaviour of A(k, cj) at kp, it can be written as, 

Mi(k) = Z^kQ{-S.k) + smooth part . (128) 

Since the slope dMiCk^/dk has a discontinuity of Zvp at kp, Paramckanti, 
et al. extracted the nodal Fermi velocity vp as shown in figure HUa). The 
doping dependence of vp together with its bare value Vp are shown in figure 
HD^b). Wc see that Fermi velocity is only weakly doping dependent, a result 
which is consistent with the ARPES data. However, this estimate of vp is 
rather inaccurate compared to the direct evaluation [146| from the quasiparticlc 
excitation energies, which will be discussed in section [T] 

Ground state expectation values also provide important information about 
the optical conductivity in the Hubbard and the t-J model. The total optical 
spectral weight -Dtot(a;) can be calculated by [11 7j 

poo 

/ duj n<j{uj) = 7r^m-i(k)n(k) = 7rAot/2 , (129) 

where m~^(k) = (9^e(k)/9k2;9kj.) is the non-interacting mass tensor. e(k) is 
the non- interacting dispersion (we set fi = c = e = 1). Since the integral in 
(|129p goes from to -l-oo, it includes contributions from the upper Hubbard 
band and is thus finite even at a; = as shown in figure HHT a). 

Paramekanti, et al. |117[ [21] emphasised that the low frequency optical 
weight, or Drude weight [202| . 

D,o^ = d^{HA)/dA\ (130) 

is more interesting, because the upper cutoff is chosen smaller than U and thus 
excludes the upper Hubbard band. In (|130p . A is the electron- magnetic vector 
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Figure 42: (a) The first moment Mi(k) of the spectral function along the zone 
diagonal, with smooth fits for k < kp and k > kp, showing a discontinuity 
of Zvp in its slope at kp. (b) Doping dependence of the nodal quasiparticle 
velocity obtained from the slope discontinuity of A/i(k). Error bars come from 
fits to Mi(k) and errors in Z. Also shown are the bare nodal velocity Vp, the 
slave boson mean field v^p{x) (dashed line), and the ARPES estimate ■yj^'^^'-' 
mm]. Model parameters: U = 12t, t' = -t/A. From [H]. 
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Figure 43: (a) Doping dependence of the total (-Dtot) and low energy (Aow) 
optical spectral weights (b) The optical spectral weight I?iow versus the nodal 
quasiparticlc weight Z. Model parameters: U = 12t, t' = —t/A. From [21] . 
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potential, which is introduced into the Hamiltonian ([T]) in terms of a Peierls 
substitution |202| . 

4aCja cl^Cja exp {i eA (R, - R^)) , (131) 

where we used A = {A, 0) and set h ^ c ^ I for simphcity. As shown in figure 
I43f a). the Drude weight Z^iow vanishes hnearly for x — > 0. This demonstrates 
that the Gutzwiller projected superconductor indeed describes an insulator in 
the half-filled limit, which can be argued to be a general property of projected 
states |117] . The VMC results for the Drude weight Uiow resemble experimental 
data in magnitude as well as in the doping dependence quite well [1171 1203] . By 
plotting I?iow versus Z (from the nodal point) Paramekanti, et al. also illustrated 
that flow oc Z, see figure HSTb). 

The Drude weight -Diow also provides an upper bound to the superfluid 
stiffness Dg, i.e. Dg < -Diow |204j . It follows that Dg ^ as x ^ in 
agreement with experiments [205] . Since the penetration depth Al is related to 
Dg by A^^ = Ane^Dg/h^c^dc , where dc is the mean-interlayer spacing along 
the c axis in a layered compound, Paramekanti, et al. |117| could also estimate 
a lower bound for Al which is again consistent with experimental data. 

These VMC calculations based on a Gutzwiller projected superconducting 
ground state describe several key features of HTSC remarkable well. The results 
are in general agreement with RMFT and confirm the usefulness of projected 
wave functions in the context of HTSC. Although restricted to T = 0, the 
above ansatz gives us some hints about the finite temperature regime. The 
superconducting order parameter $ resembles the doping dependence of Tc and 
vanishes at half-filling, while the superconducting gap (expected to scale with 
A) remains finite. This suggests that the underdoped regime exhibits strong 
pairing and that the superconducting transition may be determined by the onset 
of phase coherence (rather than the vanishing of pairing amplitude). 

6 Quasiparticle states within renormalised mean 
field theory 

Extending the RMFT to excited states requires the consideration of Gutzwiller- 
Bogoliubov quasiparticles within the t-J and the Hubbard model. These Gutz- 
willer-Bogoliubov excitations then allow for a systematic analysis of the single 
particle spectral function and explain momentum- and doping-dependent fea- 
tures in ARPES and STM experiments. Apart from these key results, in this 
section, we also discuss the renormalisation of the current carried by Gutzwiller- 
Bogoliubov quasiparticles and the consequences for the suppression of the su- 
perfluid density. We also discuss the discrepancies between different approaches 
in determining the underlying Fermi surface of a Gutzwiller projected super- 
conductor. 
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6.1 Coherent and incoherent spectral weight 

To model the spectral features of HTSC, we need to study the excited states 
of a Gutzwiller projected superconductor. In this subsection, we consider the 
transfer of spectral weight from coherent quasiparticles (QPs) to an incoherent 
background. Stimulated by STM, which reveals a striking particle-hole asym- 
metry in the spectra of underdoped HTSC [501 [53l [54] , this problem has received 
much attention recently and investigated using both RMFT [SI 11531 1206] and 
VMC methods [i EHl [13 WH [IMl [UHl [Ml [HI]. As nicely seen in the 
experiments, e.g. in figure 1(c) and 3(e) of [53] or figure [TTl the spectral weights 
for hole and particle addition show a distinct asymmetry. 

6.1.1 Sum rules for the spectral weight 

The asymmetry in the STM spectra may be explained qualitatively by consid- 
ering sum rules |206[ 12131 1214[ 12151 1216| for the one-particle spectral function, 

A(k,w) = ^(0|4jm)(m|cka|0)(5(c^ + (£„,- So)) (132a) 

771 

+ Y,i(^\cka\m){m\cljO)6iLu-iE,„-Eo)) , (132b) 

m 

with, 

/oo 
dujA(k,uj) = 1 . (133) 
-OO 

In p32a|) and (fT32b|) . we use the T = spectral representation of A(k, uj), where 
I to) arc the exact many-body eigenstates with energies Em. The ground state is 
given by m = 0, and u; is measured with respect to the chemical potential. We 
are now interested in the low energy spectral weight of a doped Mott insulator 
described by a Gutzwiller projected ground state, i.e. |0) ~ |^) = PgI^o)- 

When removing a hole from the ground state [as in (|132ap ] it is clear that 
no doubly occupied sites are created. Thus, the resulting state is situated in the 
so-called 'lower Hubbard band' (LHB) and involves only low energy excitations, 
i.e. < E„i ~ Eq U (excitation energies much smaller than the Hubbard 
U). Thus, on the hole side, the low energy spectral weight corresponding to 
momentum k and spin a is given by, 

/ dLuA{k,u;) = (0|4^Cka|0) = (nk.)vt . (134) 
<y —oc 

By summing over all spin and momenta, we obtain the total low energy spectral 
weight for the hole side, 

1 f° 1 

dujA{k,uj) = -Y,M^^n. (135) 

k,o- -^^"^ k,cr 

We note that similar sum rules can be derived for the dynamical conductivity, 
viz. the /-sum rule |217| . 
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The situation is different when adding an electron to the ground state [as in 
(|132bp ]. In such a process a part of the resulting state is located in the 'upper 
Hubbard band' (UHB), i.e. a doubly occupied site may be created. Therefore, 
we have to choose an upper cutoff ^II (located between LHB and UHB) to 
extract the low energy spectral weight. By integrating A(k, w) from to fii, we 
restrict ourselves solely to the Gutzwiller projected eigenstates out of all \m), 
and we obtairF^. 

/ duJA{k,u;)^{0\cl,aPGciJO)^{PGC^.aPGcl,PG)^o ■ (136) 
Jo 

Summing again over all spin and momenta and making use of Fourier trans- 
formation, we find the total low energy spectral weight for the electron side to 
be, 

dcMKuj) = |E(^GCk.PGcL^^G)*o (137a) 



k,o- k 



Ty.{PGCiM-ni^^)clPG)^, (137b) 



L ^ 

La 

= ^^((1 -n,_„)(l -nz,))* (137c) 

La 

= 2 • (1 - n) , (137d) 

where we used, Pg Pg{^ — ni-a) cj^ Pg (for a site /), to get (jl37bp . 

From (|135p and (|137p . we find that it is more difficult to add an electron 
to the LHB than to extract one in a doped Mott insulator. This asymmetry 
increases as one approaches half-filling. For a hole density, x = 1 — n, the total 
spectral weight one the particle side is reduced to 2x = 2(1 — n), while the hole 
side of the spectral weight is not much affected. Note that these sum rules while 
explaining the particle-hole asymmetry of the total spectral weight, tell us very 
little about the energy distribution of spectral weight within the LHB. 

Wc further note that the total spin-integrated spectral weight is 2, and the 
integrated spectral weight of the upper Hubbard band is consequently, 2 — 
n — 2 • (1 — n) = n, which agrees with the Hubbard-I approximation for the 
paramagnetic case |218j . 



6.1.2 Definition of coherent quasiparticle excitations 

To explain the distribution of spectral weight at low energies, we approximate 
the eigenstates |m) by the Gutzwiller-BogoHubov quasiparticles, equation (|103p . 
derived from RMFT [71 1206] . We formulate particle- like Gutzwiller-Bogoliubov 
QPs by, 

K+'>=^n+iPg7LI*o) , (138) 

For a more detailed reasoning leading to this step, we refer to [2061 . 
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as well as hole-like Gutzwiller-Bogoliubov QPs with the same momentum and 
spin by, 

=^n-iPg7LI*o) . (139) 

In the following, we fix the particle number N by the operator Pn and thus the 
ground state is l*'^) = Pn-PgI^'o)- To avoid confusion, we include an index N 
for the particle number in the wave function. At the level of mean field theory, 
the energies corresponding to the states (|138p and (|139[) are given by the RMFT 
excitations £'k, as discussed in section 

Using (fB5)) and ([1^ in (fH^ and (fT32b)) yields, 

A(k, uj) = Z+ul 6iu; - Sk) + Z^^vl 5{ij + Pk) + ^""(k, lj) , (140) 

with the QP weights Z^^ given by, 



and 



^k,T - ^'k - , .,N-1 |,T,N-1 \/,T,N|,T,N\ ■ 



Here, we distinguish between the QP weight Z^^ mostly used in VMC cal- 
culations and the QP weight renormalisation Z^^ often given within RMFT 
studies. In p40p . the Gutzwiller-Bogoliubov QPs lead to i5-peaked excitations, 
that are associated with the coherent peaks, e.g. seen in ARPES. For projected 
wave functions, the weight of these coherent excitations is renormalised (due to 
Gutzwiller projection) by a factor Z^. Thus, by construction, the full spectral 
weight is not exhausted by the Gutzwiller-Bogoliubov excitations, (jl38p and 
(|139p . demanding the presence of an incoherent background A"^'^(k, w). 

It is not settled yet if the asymmetry in the HTSC comes from the incoherent 
part as dictated by the spectral sum rules, or if such an asymmetry is present 
in the coherent QP spectrum [SI 11531 12061 12091 1210] . As we show below, recent 
works based on the GA support the former proposal [7l l206j . i.e. a particle-hole 
symmetric quasiparticle weight renormalisation, 

Zi. = Z+ = Z^ . (143) 

However, recent VMC calculations |209[ I211| (discussed in section 17. ip claim 
that this symmetry is exactly fulfilled only for k at the (underlying) Fermi 
surface. Therefore, zero (or very low) energy excitations would still exhibit 
particle-hole symmetry, whereas coherent excitations at higher energies could 
lead to an asymmetry in spectral weight. This asymmetry in VMC results is 
most pronounced in the underdoped region and disappears in the limit of zero 
pairing (projected Fermi sea) |210[ I211j . 
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6.1.3 Incoherent background of the spectral weight 

Next, we discuss the incoherent background of the hole spectrum. By using the 
spectral representation, (|132ap and (jl32bp . together with A(k, w) from p40p . 
we find the relation, 

(nk.>* = Z^ul + nZ , (144) 

with 



/ dujA''"'{k,uj) . (145) 
<y —oo 



Thus, the momentum distribution, {rikcr), provides the total spectral weight 
with momentum k and spin a at the hole side, i.e. the coherent weight Z^u^ 
overlaid by the incoherent background n^^^. We will calculate these quantities 
in section by the GA and show their behaviour in the first Brillouin zone. 



6.1.4 Divergent k-dependent self-energy 

In section 15. 3[ we discussed VMC results for the QP weight renormalisation at 
the nodal point k^?. These calculations show that Z ^ for x ^ 0, where 
Z = Z^ = . Before extending our considerations to all k-points, let us 
discuss some consequences for the self-energy in the half- filled limit. Due to the 
vanishing gap along the nodal direction, (0,0) (tt, tt), we can approximate 
the Green's function in the vicinity ofki?byG^^(k,a;) = w — e(k) — /x— I](k,w), 
where E = S' -I- lE". Standard arguments then lead to the results [219j . 

where the right hand side is evaluated at the node (ki?,a; = 0). Since Z — > 
for X ^ 0, |(9E'/9a;| diverges like 1/x in this limit. Due to the finite Fermi 
velocity vp (see section and section [53]) . a compensating divergence in the 
k-dependence of the self-energy with 

automatically shows up. This limiting behaviour of and Z is also experi- 
mentally observed and transcends conventional Landau-Fermi liquid behaviour, 
where the k-dependence of the self-energy is usually small. 

Equation (|147p constitutes a key experimental result for the HTSC, since 
ARPES shows unambiguously that vp const and Z ^ for x ^ [23l 
El Ell [39]. The fact that (fliT]) arises naturally within the Gutzwiller-RVB 
framework provides a strong argument for the basic premise of the theory. It 
is a consequence of the vanishing of the number of free charge carriers ^ I — n 
due to the projection close to half-filling. The number of charge carriers is, 
in contrast, ~ n and not singular within normal Fermi liquid theory. These 
considerations lead to further consequences for higher-energy features of the 
one-particle self-energy, which have been explored by Randeria, et al. |219| . 
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6.2 Calculation of the quasiparticle weight within RMFT 

To evaluate the QP weight in plT|) and pl^ within RMFT, one can follow 
[B] and use the GA for partially projected states as presented in section 13.31 
Here, we briefly discuss the mean steps required for this calculation, and refer 
to Fukushima, et al. [B] for more details. For simplicity, one may work with a 
particle excitation, 

l*k;')=^N+iPGcLl*o) , (148) 

and a hole excitation, 



l*k."'> = PN-lPGC_k-.|*0> . (149) 

Note that this redefinition docs not effect the final results since all calculations 
include norms and, 7^0-1^0) ~ cj^crl^o) ^ c_k-cr|^o)j for a BCS wave function 
l*o). 

The first step in the calculation of the QP weight is the determination of the 
norms 

of the excitations 1^'^^^)- Invoking GA for partially projected states (see section 
13. 3p . one finds [B], 

where gt = (1 — 7i)/(l — n^-), = and Ti^cr{^t.a^ka)^o the momen- 

tum distribution function in the unprojcctcd wave function. It should be noted 
that the above result was derived for the non-magnetic case tt-o- = n-f = = 
n/2. 

After calculating the normalisation, one can use the same techniques to get 

0, 



(v1/o|ci,,PgPn+icL^nPg|^o) 

(M'olcL^N-lPGCk.i'GPNl*^) 



N, 



ff,(l-<J, (151) 



Ka , (152) 



G 



for the numerators in the equations for the QP weights, p4ip and (|142p . Using 
(|15ip . (|152p . and the normahsations in (|14ip and (|142p . we find the the QP 
weights of particle- and the hole- like excitations, 



K< = .JmJl.^x,.::..... ^ g^(i-<j , (153) 
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Figure 44: Rcnormalisation of the Gutzwiller-Bogoliubov nodal quasiparti- 
cle as a function of doping x. The model parameters are t = — 1'/4 and U = 12t. 
RMFT results for the Hubbard and the t-J model are compared with VMC data 
for the Hubbard model (from [21]) and with the SBMFT result in the t-J model. 
From [7]. 

respectively. Since n^^ = ^'k ^ ^ ^ '^k' follows w w gt, which vanishes 
at half filling n ^ 1. 

The above results show that within the Gutzwiller approximation, the co- 
herent QP weight does not cause a particle- hole asymmetry, i.e. Z^ « Z^ . 
It seems therefore that the asymmetric DOS observed in STM can only be ex- 
plained by the incoherent spectrum of Gutzwiller projected superconductors. A 
symmetric spectral weight for coherent QP excitations is also obtained in calcu- 
lations for the Hubbard model (include transformation e^*'^) [71 1206] . However, 
the RMFT results for Z^ = Z^u^ and Z^ = Z^v^. do not exactly match recent 
VMC calculations [2091 [2101 EIH [212] , which directly evaluate Z^ and Z^ (see 
section FTT]) . Nevertheless, the general doping dependence of above QP weight 
qualitatively agrees with VMC results and with the coherent weight seen in 
ARPES measurements [7]. 

6.3 Quasiparticle weight for the Hubbard model in the 
strong coupling limit 

In the previous subsection, we illustrated how one can determine the QP within 
the GA. Here, we follow [T] 1206] and extend this calculation to the Hubbard 
Hamiltonian, in analogy to the extensions of the RMFT discussed in section 
14.31 By using a re-transformed ground state, |^') = e^^^ PgPnY^i^) as well as 
re-transformed excited states. 




(155) 
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Figure 45: (a) Quasiparticle weight Z\^n^^ and (b) momentum distribution 
(rikcr) of the Gutzwillcr-Bogohubov quasiparticle for different doping x] The 
corresponding Fermi surface, = 0, is shown in the inset of (a). The model 
parameters are t ~ — 1'/4 and U = 12t. From [7]. 



we can systematically study the QP weight renormalisation within the Hubbard 
model in the strong coupling limit. Evaluating the canonical transformation 
g-iS order 0{t/U) gives the following particle-hole symmetric QP weight 
renormalisation [3j -^k = = , 

withe^ = 2t(cosfca;-|-cosfcj,) + 4t'cosfc^cosfcj,. Equation (|156p also includes cor- 
rections from the next nearest neighbour hopping term t' . The renormalisation 
Zk of the nodal QP weight is plotted as a solid line in figure [44l and agrees well 
with VMC results for the Hubbard model |1T]. The dashed line corresponds to 
the RMFT result for the t-J model, Zk = gt, which is compared to the dotted 
line, Zk = X, from slave boson mean field theory (SBMFT). 

The spectral weight of the coherent peak, measured in ARPES, is related to 
the QP weight Z^ = Z^n'^^^; it is shown in figure [45Fa) along the directions, 
(0,0) (7r,0), (tt, 0) (7r,7r), and (tt, tt) (0,0) for different x. As seen in 
the figure, the QP spectral weight is severely modified by Gutzwiller projection. 
It decreases with doping, and vanishes at half filling. This causes a shift of spec- 
tral weight to an incoherent background as seen in the momentum distribution 
function, (rikcr) ~ Zkt'l + n^^^ + 0{t/U)'^. While the first term corresponds 
to the coherent QP weight, the second gives the distribution of the incoherent 
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Figure 46: Doping dependence at the antinodal point, k = (tt, 0): (a) QP renor- 
malisation the unrenormaliscd QP weight, n^^ = v'^, and the renormahsed 
coherent QP weight Z\^i)^] (b) coherent weight Z\^v^ compared with the ex- 
perimentally determined Superconducting Peak Ratio (SPR, ratio of coherent 
spectral weight relative to the total spectral weight) for Bi2212 (37]; see also 
figure [51 The model parameters are t = —t' /A and U ~ 12t. From [7]. 



part. One obtains [7], 
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which is a smooth function of k, where At- and are the pairing and hop- 
ping amplitudes between nearest and next nearest neighbour sites, r = (±1,0), 
(0,±1), (±1,±1), as defined in section ID Results are shown in figure llSl b) 
. The incoherent weight is spread over the entire Brillouin zone, and overlies 
the coherent part from the Gutzwiller-Bogoliubov quasiparticles. At half-filling, 
all weight becomes incoherent. These results are in qualitative agreement with 
calculations for the t-J model (recovered by neglecting the </J7-corrections in 
above equations). 



6.3.1 Non monotonic behaviour of the QP weight at (tt, 0) 

Here, we consider the coherent QP weight ^k^^ at the antinodal point k = (tt, 0) 
within the Hubbard model in the strong coupling limit {U = 12t). The RMFT 
theory predicts a non monotonic behaviour as a function of doping, shown in 
figure [45F a) and figure [46l This effect arises from a combination of the effects 
due to the Gutzwiller projection and to the topology change [see insert of figure 
l45l a)] of the underlying Fermi surface (FS); figure llGT a) illustrates this clearly. 
While the QP weight renormalisation, Zk, increases with increasing doping, 
= '^k' decreases due to the topology change, which occurs at x « 0.15 — 0.20 
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Figure 47: The magnitude of the current |j(k)| along the nodal direction, (0,0) 
to (tt, tt). (a) VMC calculations on a 10 x 10-lattice (data taken from [208j ) and 
(b) RMFT results are compared for the doping levels x = 0.01,0.05,0.09,0.17 
(increasing magnitude). From [TO] , 

for our choice of hopping parameters [i! = —tj^. The change of the FS seems to 
be a generic feature of hole doped Cuprates [22011221] , although the exact doping 
concentration x, for which this occurs, is sensitive to the ratio between various 
hopping parameters. The combined effect of strong correlations and topology 
change leads to a maximum of the QP weight for the doping level, cc, at which 
the underlying FS changes topology. Indications for such a behaviour have been 
seen in ARPES [37l [38] . Feng, et al. [37] extracted the superconducting peak 
ratio [SPR, illustrated in figurcHSTb)] which is proportional to the coherent QP 
spectral weight, v^. They found that the SPR increases with small a;, attains 
a maximum value around x w 0.2 where it begins to decrease. Ding, et al. [55] . 
reported similar results from ARPES. In figure I^CT b). the SPR experimentally 
drops below the theoretical prediction for underdoped samples. This is likely to 
be the effect of inhomogeneities and of the resulting gap variations [54] , which 
cause a strong scattering of quasiparticles near the antinodes. 

Although the topology change does not influence the stability of the su- 
perconducting state within RMFT, the superconducting pairing parameter $ 
(related to Tc) and the QP weight Zi^^v^ show some similarity as a function of 
doping. However, we emphasise that this similarity does not result from any 
direct relation between these two quantities. 

6.4 Quasiparticle current renormalisation 

An important issue in the theory of the high temperature superconductors are 
the properties of the nodal quasiparticle (NQP) excitations, in particular the 
renormalisation of the respective quasiparticle current [222] and their role in 
suppressing the superfluid density ps- As pointed out by several authors [223[ 
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Figure 48: Magnitude of the nodal current as a function of doping. RMFT 
results arc compared with the VMC calculations on a 20 x 20-lattice (data 
taken from |208| ). Insert: RMFT result for the current rcnormalisation factor, 
oc = |jnodai|/i'F, as a function of doping x. From [10]. 



I224[ I225j , the proliferation of NQPs at finite temperatures decreases ps [T) |223| , 

^ = ^ - f T , (158) 

where vp and V2 are the NQP velocities in the longitudinal and transverse 
directions respectively, and the zero temperature superfluid density. The 
renormalisation factor a (also called effective charge |225j ) relates the current 
carried by the quasiparticle to its velocity, 

j(k) = -eav(k) . 

Assuming that superconductivity is destroyed by thermal NQPs , Tc is deter- 
mined by simply setting (jl58[) to zero, i.e. determining the temperature at which 
the superfluid density vanishes [2231 1224j . The behaviour of Tc as a function of 
doping is then governed by the doping dependencies of the various quantities in 
(|158p . The latter can be calculated within the framework of the RVB theory. 
Numerical [5T] calculations show that — > as x — > 0. The nodal velocity 
vp is approximately constant [7], whereas the transverse velocity V2 increases 
as the insulator {x = 0) is approached. The situation is rather unclear for the 
renormalisation factor a. While some theories argue for a constant a [225], 
recent experimental (measurement of the superfluid density [226]) as well as 
theoretical results [TO] I208j seem to support the conclusion that a decreases as 
X ^ 0. 
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To clarify this issue, Edegger, et al. [TU] used RMFT to calculate the current 
renormalisation for the t-J model with J = i/3. For the supcrfluid density at 
zero temperature, RMFT yields a doping dependence of, 



where we used |202j . 



pio)^g,^J^, (159) 



'^U{cl+r,aC^,a+claC^+r,a)) , (160) 



and evaluated (|160p . invoking the Gutzwiller approximation. Here, we used 
T = x,y and neglected corrections due to the re-transformation e"''^ of the 
wave function to the Hubbard model, i.e. we set e~''^ = 1. Using linear re- 
sponse theory for the supcrfluid density |202j and restricting ourselves to low 
temperatures, we recover ()158p within RMFT [10] . The renormalisation factor 
a can be derived by considering the current carried by the Gutzwiller projected 
Bogoliubov quasiparticlc states l^'ka), 

j(k) =ie{Y^ Uj (^4^ Cj-^ - c]^^ Ci^a^)'ii,^ ■ (161) 

By invoking the Gutzwiller renormalisation scheme, we find, 

j(k) = -e5,^e°(k), (162) 

where e''(k) is the unrenormalised tight binding dispersion relation; again we 
set e^''^ ~ 1 for simplicity, i.e. we neglect any </J7-corrections in (jl61[) and 
l|162|) . Combining (|159p and (|162[) allows us to extract a. At the nodal point, 
one finds a = gtVp/vp, where Vp is the unrenormalised Fermi velocity. The 
results are shown in figure ST] and figure [35] along with VMC data taken from 
|208j . As can be see, both methods are in excellent agreement and show that 
the renormalisation factor a ^ 0, as a; ^ 0. Since the x dependence of the 
supcrfluid stiffness can be obtained experimentally, it is important to study 
dps{T)/dT (X o?VFjv2. We show the results for this quantity in figure HWa) . 
Note that vf/v2 oc vp/Asc already shows a significant x-dependence and may 
explain the experimentafly observed doping dependence of dps{T)/dT [226] . 
However, multiplication by leads to a slope dps{T)/dT that vanishes as 
X ^ 0, i.e. as X 0, the effective NQP charge vanishes faster than the 
supcrfluid density does. Therefore, we get a meaningless estimate for by 
setting (|158p to zero as shown in figure 

This problem was noted by Lee and Wen [2231 1224] in the context of the 
U{1) gauge theory of the t-J model. They argued that an SU{2) formulation 
may resolve the problem, yielding a constant a. However, a constant a does 
not completely agree with the experimentally observed x-dependence of the su- 
perconducting dome either [maximal Tc at x w 0.08, see figure HHT b)]. There 
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Figure 49: (a) Doping dependence of vp/v2 and o^vfIvi from RMFT with 
t' = — 0.2t and U = Ylt in the large U Hubbard model, (b) Doping dependence 
of Tc from setting pSSp to zero: (i) a[x) from RMFT (ii) oi[x) is set equal to 1 
by hand. From [10] . 



are several possible reasons for the discrepancy. It may be that the RMFT 
result for a is indeed correct, in which case, the issue can be resolved by more 
experiments explicitly extracting a in the underdoped regime. This would au- 
tomatically mean that is not determined by NQPs, i.e. (jl58p . and one needs 
to look for other possibilities such as vortex proliferation as mechanisms that 
set the scale for T^. 

Another possibility is that the theoretical framework behind the Gutzwillcr 
RVB theory misses a crucial ingredient in the derivation of (|158p and the calcu- 
lation of the effective current renormalisation a. Indeed, the applicability of the 
standard Kubo formula for ps 2021 in a projected Hilbcrt space may be ques- 
tioned and one needs to reexamine this calculation carefully to check whether 
(|158p is indeed correct. 

A recent more phcnomenological approach argues that the overall tempera- 
ture dependence of the superfluid density at low dopings is well described by a 
three-dimensional strongly anisotropic weakly interacting Bose gas |227| . How- 
ever, more work is necessary to connect such phenomenological models to the 
RVB theory we outlined so far. 

6.5 Determining the underlying Fermi surface of strongly 
correlated superconductors 

The underlying Fermi surface (FS) in the HTSC was studied recently by us [9] 
and Sensarma, et al. |228j . These results clarify the notion of a FS in a super- 
conducting state and what does it mean when we say that ARPES measures 
the FS of a superconductor. 

In the case of the HTSC, due to the large superconducting gap (pseudogap) 
below (above) the superconducting transition temperature, an FS can be defined 
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only along the nodal directions (the so-called Fermi arcs ^51 [i^ liT] V 
The full 'underlying FS' emerges only when the pairing interactions are turned 
off, either by a Gedanken experiment, or by raising the temperature. Its exper- 
imental determination presents a great challenge since ARPES is more accurate 
at lower temperatures. Therefore, it is of importance to know what is exactly 
measured by ARPES in a superconducting or in a pseudogap state of the HTSC. 



6.5.1 Fermi vs. Luttinger surface 

We follow [3 and begin by highlighting the differences between a Fermi and 
a Luttinger surface. The FS is determined by the poles of the one electron 
Green's function G{\i,uj), viz. by ReG{k,uj = 0) = ±oo |229| . The Luttinger 
surface is defined as the locus of points in reciprocal space, where the real part 
of the one particle Green's function changes sign |230] . In the Fermi liquid 
state of normal metals, the Luttinger surface coincides with the FS. In a Mott- 
Hubbard insulator the Green's function changes sign due to a characteristic 
1 / w-divcrgcnce of the single particle self energy [231i I232[ 1233] at momenta k of 
the non-interacting Fermi surface. In the HTSC the gapped states destroy the 
FS but only mask the Luttinger surface. Hence, it seems natural to relate the 
Luttinger surface of the superconducting and of the pseudogap states with the 
concept of an 'underlying FS', and ask if such a surface can be determined by 
ARPES. 

The single particle Green's function is given by, 

^ ' ' ^ uj - {En - Eo) + iO+ ^uj + {En - £^0 + ^0+ ' ^ ' 

where En are the eigenvalues corresponding the eigenstates jji) of the Hamilto- 
nian; the ground state and its energy are given by |0) and i^o, respectively. In 
order to perform explicit analytic calculations one can approximate the coherent 
part of (|163p by the RMFT results for the Hubbard model (sec section and 
section |673)) . In analogy to section [6.1.21 for the spectral function A(k, w), we 
can use, Z^u\ = |(ri|4^|0)p, Z^v\ = | (7i|cfc^|0) and, E^ = En - Eq. Thus, 
the RMFT result for the coherent part of the Green's function becomes 

g(k,.)^ Jfl.^, . (164) 

Within RMFT the elementary excitations in the superconducting d-wavc ground 
state are given by the dispersion relation. 



E^ = + Ai , (165) 

where and Ak arc determined by (|108p and p09p . respectively. Evaluating 
ReG(k,w = 0) by (UMI) one finds 

RcG{k,co = 0) = ^{vl-ul) = -^^^, (166) 
iik E^^ 
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Figure 50: The zero frequency spectral intensity (deduced from the inverse of 
i?k! which was determined by RMFT with model parameters t' — —t/A and 
U ~ 12t) in the first Brillouin zone for hole dopings x = 0.05 (left) and x = 0.25 
(right). The colour coding blue/red corresponds to the low/high zero frequency 
spectral intensity. The ridges of maximal intensity are indicated by the (dashed) 
red and (dashed-doted) orange lines respectively, the Luttinger surface by the 
black line. From [5]. 



where the right hand side follows from the mean field relation, v^= 1 — = 
(1 - Ck/£'k)/2 [see ^] . The poles of Re G{]i, w = 0), which determine the FS, 
are therefore given by 

E^ = Q. (167) 

However, for a d-wave superconductor, equation p67p is fulfilled only at the 
nodal points; consequently a FS is well defined solely at these points. Alterna- 
tively, one can consider the Luttinger surface, defined by sign changes in the 
Green's functions at = 0. From (|166p . sign changes are found whenever 

Ck = . (168) 

From above equations, we conclude that the Luttinger surface is determined by 
the condition = 0, which is also the definition of the normal state FS when 
Ak = 0. 



6.5.2 Fermi surface determination 

There are several ways to determine the FS in practice. However, these meth- 
ods do not coincide with the underlying FS, viz. the Luttinger surface, in the 
HTSC due the large superconducting gap (or large pseudogap) in the under- 
doped regime. 

To demonstrate this fact, we follow [9l I228j and discuss the so-called 'maxi- 
mal intensity method' in more detail. In this approach the intensity of ARPES 
spectra at zero frequency is used to map out the underlying FS. This quantity 
is determined by A(k, = 0) = — ilmG(k,cj = 0), which becomes 

(169) 



if one replaces 0+ by a finite broadening Fk in (|164p . The Fk is determined 
both by the experimental resolution and the width of the quasiparticle peak. 
When the momentum dependence of Fk is small compared to that of (as is 
usually the case), the maximal intensity is given by the set of momenta Kk. for 
which i?k is minimal. 

This method in determining the underlying FS was examined in [5] by cal- 
culating (|169p within RMFT for a strongly correlated d-wave superconducting 
state. All calculations in [9| were done with model parameters for HTSC using 
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RMFT [3 HO] , for which the quasiparticle dispersion _Ek retains the form of 
(|165p . Figure [50l shows results for the spectral intensity at zero frequency as 
well as the locus of the Luttinger surface, where the former is deduced from the 
inverse of i?k- 

For large hole doping, x = 0.25, the superconducting gap is small and the 
Luttinger surface is close to the points in momentum space for which the zero 
frequency intensity is maximal. But for smaller doping, x = 0.05, the gap is 
substantial and the Luttinger surface deviates qualitatively from the maximal 
intensity surface due to the momentum dependence of Ak (see ridges in figure 
[50]) . It follows that when the gap or the pseudogap is large, the criterion of 
maximal spectral intensity alone does not suffice to identify the correct FS and 
it is necessary to supplement the analysis of the zero frequency ARPES intensity, 
(|169p , with a dispersion relation such as (|165[) . These considerations explain why 
the (outer) maximal intensity ridges seen in ARPES (at low temperatures in 
the underdoped regime) may yield an underlying FS whose volume is too large. 
In particular, this effect is seen in Ca2-xNaa;Cu02Cl2 |234j . which also exhibits 
quite a large pseudogap |235| . 

As discussed by Gros et al. [9], even larger deviations from the underlying 
FS are present in the 'maximal gradient method'. This method is based on 
the fact that the FS is given by the set of k-values for which the momentum 
distribution function nk shows a jump discontinuity. When this discontinuity 
is smeared out, say, by thermal broadening or a small gap, the gradient of rik, 
jVk'n.kl, is assumed to be maximal at the locus of the underlying FS. However, 
this method is very sensitive to the presence of even small gaps [9] and cannot be 
used to determine the underlying FS unambiguously from numerical [236|, I237j 
or ARPES data [238|,I239] . Furthermore, we note that even the Luttinger surface 
in the HTSC can slightly violate the Luttinger count. This surprising result is 
discussed in [9] and |228| in more detail. 

Bieri and Ivanov [21 Ij recently proposed an alternative definition of the 
underlying Fermi surface kp- by the condition Z^^ = Z^^, viz. that the quasi- 
particle and the quasi- hole weight coincide at the FS, as they do for a Fermi 
liquid state (see section [01 for the definition of Z^). This definition also agrees 
with the Luttinger surface, ^k = 0, within RMFT. However, when considering 
a Gutzwiller-projected superconducting state within the VMC technique devia- 
tions from the Luttinger surface are observed |211j . This deviations stem from 
the asymmetry between Z^ and Z^ , which shows up in the VMC calculations 
only. 

6.5.3 Renormalisation of the Fermi surface tov^rards perfect nesting 

The presence of strong electron-electron interactions also changes the geometry 
of the Luttinger surface close to half filling. The Cu-0 planes of the HTSC 
are characterised by a nearest neighbour (NN) hopping parameter t « 300 meV 
and a next nearest neighbour (NNN) hopping parameter t' ~ —t/A. These 
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hole doping x hole doping 



Figure 51: (a) Renormalisation of the next nearest neighbour hopping ampH- 
tude, t' t', as a function of hole doping concentration x for various values 
of bare t'. All effective i' are renormalised to zero at half filling by the large 
Coulomb repulsion. We highlight the region for which we expect the supercon- 
ducting d-wavc state to become unstable against antifcrromagnctism (AFM) 
due to the nearly perfect nesting of the Luttingcr surface. RMFT calculations 
with U = 12t. From (b) Determination of t'^^ by a fit to the maximal inten- 
sity surface, see figure [501 Wc compare RMFT calculations with Cluster-DMFT 
(C-DMFT) calculations from CiveUi, et al. [240] for t' = -0.3t (indicated by 
the black line) and U = 16t. 



parameters are the bare parameters, and determine the dispersion relation, 

ek = —2t{ cos kx + cos ky ) — 2t' ^ cos{kx + ky) -f cos,{kx ~ ky) ^ , (170) 

in the absence of any electron-electron interaction. On the other hand, true 
hopping processes are influenced by the Coulomb interaction (here U = \2t) 
leading to a renormalisation of the effective hopping matrix elements, 

t i^i{U), t' i'^i'iU). (171) 

One can extract t and t' from the RMFT dispersion in (|108p . see section l473l 
and finds close to half-filling f cx J = At"^ /U and i' 0, i.e. the NNN hopping 
is renormalised to zero. This behaviour is illustrated in figure [STT a). The 
resulting Luttingcr surface renormalises to perfect nesting. A similar behaviour 
has been observed in recent variational studies of organic charge transfer-salt 
superconductors [94] . 

At half filling the Hubbard model reduces to a spin-model with NN J = 
4t'^/U and a frustrating NNN J' = Mt')'^/U. The ground state wave function 
obeys the so-called Marshall sign rulcoin the absence of frustration, J' = 0, viz. 
when the underlying Fermi surface is perfectly nested by the reciprocal magnetic 

^® Marshall |241) showed that the ground state of the spin-^ Heisenberg Hamiltonian on 
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ordering vector Q = (tt, tt) (in units of the inverse lattice constant). Hence, any 
deviation from the Marshall sign rule as a function of the frustrating J' can be 
used to determine the degree of effective frustration present in the ground state. 
We emphasise this is a qualitative statement of the ground state wave function. 
A numerical study has found, that the Marshall sign rule remains valid even for 
small but finite J', viz. the effective frustration renormalises to zero |242| . Such 
a behaviour is in agreement with the results presented in figure [STf a). 

However, we note that the renormalisation to perfect nesting was not seen 
in Cluster-DMFT (C-DMFT) studies, e.g. by CiveUi, et al. [HD]. We beheve 
that these discrepancies stem from the way of fitting the effective NNN hopping 
t'^g. Within the C-DMFT study of Civelh, et al. the FS is determined very 
similarly than within the maximal intensity method. Therefore, deviations from 
the effective t' of the Luttinger surface [shown in figure ISlT a)] are unsurprising. 
In figure ISlTb) , we compare the effective NNN hopping t'^g determined by a 
fit to the maximal intensity surface (see figure EQI outer ridges) with the C- 
DMFT results [240] . Both methods show a qualitatively very similar doping 
dependence. When approaching half-filling \t'^g\ first decreases, but then starts 
to grow rapidly in the undcrdopcd regime. We associate this effect with the 
increasing influence of the d-wave gap in the maximal intensity surface at small 
doping. Above considerations show that the determination of the underlying 
FS in the HTSC is a tricky task, where special care is required when comparing 
data from different approaches. 

7 Quasiparticle states within the variational Monte 
Carlo scheme 

VMC calculations for the QP weight in the t-J model only agree qualitatively 
with the approximative RMFT results. Minor deviations from the RMFT stud- 
ies may explain a contribution of the coherent excitations to the distinct particle- 
hole asymmetry seen in the STM spectra. Apart from the QP weight, we also 
discuss excitation energies determined by VMC calculations, which match well 
with previous RMFT results. 

7.1 Direct calculation of the quasiparticle weight 

RMFT together with GA is an useful tool to analyse QP features in strongly 
correlated superconducting states. However, the RMFT and GA are approxi- 
mate methods and it is desirable to check their predictions numerically by VMC 
calculations. This consideration motivated several authors |207|, 12081 12091 12101 
Uni Hm to calculate the QP weight, piT|) and directly by evaluating 

appropriate expectation values within the projected wave function j^"). These 

any bipartite lattice will be a singlet. Furthermore, the ground state wave function picks up 
a sign whenever two antiparallel spins from different sublattices are interchanged. This is the 
Marshall sign rule. 
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Figure 52: VMC result for Gutzwiller projected c?-wave BCS state on a 18 x 18 
lattice with 42 holes (x ~ 0.13) and A/t = 0.1. (a) Momentum distribution 
function (rik). (b) Total QP weight Zl°\ (c) QP weight in the (0,0) - (7r,7r) 
direction, (d) QP weight in the (0,0) - (0, tt) direction (total add Z^, 

remove Z^). Results correspond to the t-J model since the re-transformation 
of the wave function was neglected. From |209j . 



VMC studies confirm the RMFT prediction, that the QP weight decreases to- 
wards half-filling, where it finally vanishes. However, as we will show below, the 
VMC results reveal some limitations of the RMFT concerning the determination 
of the detailed doping- and k-dependence of the QP weight. We note that most 
of the VMC calculations presented below do not include a re-transformed trial 
wave function and describe observables in the t-J model. These calculations 
can be directly compared to the RMFT results from section [^21 

To calculate the QP weight within the VMC scheme, most authors use two 
helpful exact relations for Gutzwiller projected wave functions. First, one finds 
for the QP weight Z^^ of electron-hkc excitations, that [207 1 [208 1 [209 1 [2T0 l [2TT] . 

Zi,^^^-{n^.h^, (172) 

can be derived without any approximation and assumption. Thus Z^^ can be 
calculated from the momentum distribution of the ground state j^"^) [2071 1208] . 

For the QP weight Z^^, there is no exact relation corresponding to (|172p . 
However, several authors showed [2101 12111 1212| that Z^^ and Z^^ combined 
satisfy the exact relation, 



-k-o- ^kcr ^^N I ,J,N^^,5N-2 I ^N-2|) 



^^k-o • ^kff - /,T,N I ,T,N\ /,T,N^-2 I ,T,N^2'\ = . (173) 



This relation is very useful, because the matrix elements contributing to 
only involve ground states with different particle numbers. The quantity 
is closely related to the off-diagonal long-range order in the pairing correlation 
and can be calculated in a straightforward way by VMC techniques (see, e.g. 
|210j ). Equation (|173p was also confirmed numerically |212j . However, we note 
that (|172p and (|173p are only valid for the projected wave functions Pg\'^), 
and cannot be used for the re-transformed wave function e~'^^ Pg\'^) , since the 
canonical transformation e~*^ does not commute with the electron and the 
projection operators. 

7.1.1 Momentum dependence of the quasiparticle weight 

VMC results for the QP weights Z^ (adding an electron) and Zj^ (removing 
an electron), together with the total weight, = Z^ + Zj^, are summarised 
in figure [52l These calculations show that is continuous over the whole 
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Figure 53: The momentum distribution (rik) along the nodal direction k = {k,k) 
for the Hubbard [black squares in (b)] and the t-J model [white squares in 
(b)] from (a) RMFT and (b) VMC calculations, respectively. The calculations 
are based on the full t-J Hamiltonian ^ with t' = —t/i and U = 12t at a 
doping level x = 0.05. Expectation values for the Hubbard model are evaluated 
within a re-transformed wave function, see (jllObp . whereas these corrections are 
neglected in the t-J model. The RMFT calculations are based on the results 
from section 16.21 and section 16.31 the VMC data are taken from [117] . 



Brillouin zone, thus supporting the idea that = at the (underlying) Fermi 
surface |209j . However, away from the Fermi surface, figure[52]also exhibits some 
deviations from the simple RMFT calculations [Z^ = gtu^ and Z^ = gtv^ with 
gt = 2.t/(1 + x)]. For instance, inside the Brillouin zone and along the nodal 
direction, RMFT gives a constant QP weight Z^ (since (nko-) — is constant 
along the nodal direction in the t-J model, see figure [55]) whereas the VMC 
calculations (green triangles in figure \^ c] . see also |211j ) clearly show a non- 
constant behaviour. 

In the absence of a superconducting gap the quasiparticlc weight at the Fermi 
surface is determined by the jump in the moment distribution (jiko-) j as discussed 
in section [5.3.2l Furthermore, Z^ is generally related via (|172p to (riko-) for the 
t-J model. Due to this relation between (riko-) and Z^ we rc-considcr the mo- 
ment dependence of (riko-) in the VMC and the RMFT calculations. In figure 
[53] we show RMFT as well as VMC results for the moment dependence of (nua) 
along the nodal direction determined within the Hubbard and the t-J model, re- 
spectively. We note that expectation values for the Hubbard model are obtained 
by applying a re-transformed wave function e^''^|^'), which can be evaluated in 
order 0{t/U); see (|110bp in section H31 On the contrary, the re-transformation 
e~*'^ is neglected for the calculation of observables in the t-J model. Figure 
[53lal shows that, except for the jump at the Fermi point k^?, the RMFT gives 
a constant (nkcr) along the nodal direction for the t-J model. However, VMC 
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Figure 54: The doping dependence of average QP weights Z~^^ for removing 
an electron in a d-wave state obtained by VMC calculations (12 x 12 lattice, 
t' = 0) and by RMFT, respectively. The results are for the t-J model (no 
re-transformation of the wave function). The squares (triangles) are the VMC 
resuhs for Z'^ {n'lH = riavo - Z~^) with riavc = l/iEk("-ka)* = (1 - x)/2. 
The dashed and dotted lines without data points represent results by RMFT. 
From piU] . 



calculations at the same doping level {x = 0.05) and for the same model parame- 
ters exhibit a non-monotonic behaviour near the Fermi point, see white squares 
in figure [55f b). This eS'ect comes from the correlated hopping nature of the 
electron in the projected Hilbert space and is not obtained within RMFT. This 
result also explains the origin of the discrepancies between the RMFT and the 
VMC methods in determining the quasiparticle weight and reveals some limita- 
tions of the RMFT in calculating momentum dependent quantities. However, 
including the re-transformation of the wave function for the Hubbard model, 
removes the non-monotonic behaviour of (nko-) in the VMC data, figure [53{b). 
Thus, RMFT and VMC are in better qualitative agreement, when (nkcr) is cal- 
culated within the Hubbard model, compare solid line in figure ISSl a') with black 
squares in (b). 

7.1.2 Doping dependence of the mean quasiparticle weight 

Some discrepancies between VMC and RMFT in the doping dependence of the 
coherent QP weight have been discussed by Chou, et al. |210j . The authors 
calculate the average coherent QP weight for removing an electron, 

k 
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Figure 55: Dispersion E in the nodal direction for the 2D t-J model with 
J/i=0.3 and t'/t=-0.2 at a-=:0.099. (a) FuU dispersion for L=162 (triangles) 
and 242 (eircles). The electron removal (addition) spectrum is denoted by open 
(solid) symbols. The dashed lines are tight binding fits, (b) Same as (a) but 
focusing on the excitations near Ep. In addition to the data for L=162 (open 
triangles) and 242 (open squares), results for L = 1250 (solid squares and circles) 
are also plotted. From [146j . 



by the VMC scheme and compare it with the RMFT results. As shown in 
figure EH VMC calculations give a significantly larger coherent QP weight than 
RMFT at the hole side, which is directly related to a reduction of the (average) 
incoherent background n"^^ by the same amount. 

On the other hand the average QP weight for adding an electron, 

~ I 1 \ ^ ~ I 1 + X 1 \ ^ . , 1 -\- X 1 X . . 

Kvc = j2^^Wa^ — J 2^("k<T)* = Y~ ^ ^ ' 

k k 

is exactly the same in the RMFT and the VMC scheme, where we used (|172p in 
(|175p . Thus, it was argued |210j . that the increased coherent weight at the hole 
side seen in the VMC calculations, can explain the particle- hole asymmetry 
in the tunnelling experiments. However, considering the large asymmetry in 
the experiments and the predictions from sum rules, it is likely that at least at 
higher energies a considerable part of the asymmetry is caused by the incoherent 
background. 

7.2 VMC calculations for the quasiparticle energy 

In the previous section, we discussed how the spectral weight of Gutzwiller- 
Bogoliubov QP excitations can be determined directly using VMC and how such 
results compare with RMFT. Now we turn to the excitation energies E\^, of the 
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QP. Here again, RMFT results can be checked by directly calculating the energy 
corresponding to the excited state ct^), equations (|138|) and (|139p . within 
the t-J model. Subtracting the ground state energy, we obtain the excitation 
energy, 

E={Ht-j)^.±i- {Ht-j)^. (176) 

k.cr 

Wc discuss now the VMC calculations of Yunoki. et al. |146j . who also in- 
cluded a Jastrow factors into the wave functions to improve the ground state 
energy. Figure [SS] illustrates a typical dispersion along the nodal direction ob- 
tained by determining ~ \E\ for every k-point separately As shown in the 
figure, a tight-binding dispersion fits well to the numerical data, and it is pos- 
sible to extract interesting quantities like the nodal Fermi velocity vp or the 
nodal Fermi point |ki?|. 

By repeating the calculation from figurc[55]for various electron densities, one 
can determine the doping dependence oi vp and |ki?|. Figure ISHl a) illustrates 
that the Fermi velocity only slightly decreases when approaching half-filling as 
already seen from RMFT ([7], section [O]). The results of Yunoki, et al. also 
agree with previous VMC calculations utilising the moments of the spectral 
function ([H], section [53)1 . as well as with ARPES experiments [551 11561 11571 
115811159] . In figure lSCT c). we see the doping dependence of the nodal Fermi point 
jk^'l, which matches experimental and RMFT predictions. The renormalised 
band width W is given in figure ISSTb): it is tightly related to vp- Figure ISCTd) 
illustrates a comparison between the Fermi velocity vp and the unrenormalised 
velocity w^, revealing the strong renormalisation effects due to the Gutzwiller 
projection. However, it is important to note that, in contrast to the QP weight 
Zk, the Fermi velocity does not vanish in the half- filled limit 

While Yunoki et al. |146j only considers the nodal dispersion, Yunoki re- 
cently extended these VMC calculations to the whole Brillouin zone |212| . His 
results agree quite well with the RMFT dispersion, giving further support to 
the Gutzwillcr-Bogoliubov QP picture. To conclude this section, VMC calcu- 
lations for the spectral weight and the QP excitations are in good qualitative 
agreement with RMFT. This is important because it shows that the simple ana- 
lytical approach of RMFT together with the GA can be used reliably. Two key 
features emerge consistently from these two approaches; a finite and constant 
Fermi velocity contrasting with a vanishing QP weight in the half-filled limit. 

8 Summary and outlook 

In this review, we attempted to summarise the basic idea of using Gutzwiller 
projected wave functions in the description of high temperature superconduc- 
tivity. Projected wave functions provide a straightforward implementation of 
the RVB picture wherein superexchange leads to pair correlations and doping 
the Mott insulator leads to a superconducting ground state. 

Projected wave functions can be studied both analytically and numerically. 
A superconducting state with d- wave symmetry arises as the best variational 
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wave function within the Gutzwillcr -RVB theory. Incorporation of antifcrro- 
magnctic order and next-nearest neighbour hopping then allows for a quantita- 
tive description of the Cuprate phase diagram within the t-J and the Hubbard 
model. Sophisticated variational Monte Carlo calculations (VMC) give detailed 
information about the size of the antiferromagnetic region and the stability 
against phase separation. These ground state properties seen in the VMC tech- 
nique were recently confirmed by various quantum cluster methods, lending 
further support to the Gutzwiller-RVB picture. 

Besides the VMC technique, the effect of projection can be treated by Gutz- 
willer approximation, which then allows for a formulation of a renormalised 
mean field theory (RMFT). The RMFT results agree in general with VMC 
calculations and provide systematic analytic expressions for doping-dcpcndcnt 
features. 

Within the Gutzwiller-RVB picture, high temperature superconductors are 
viewed as doped Mott insulators, i.e. restriction to single occupied orbitals due 
to strong correlation effects. That causes a significant decrease in the mobility 
of electrons (holes) near half-filling as correctly described within above micro- 
scopic calculations. The resulting renormalisation of the kinetic energy explains 
the decrease of the superconducting order parameter, of the superfluid density, 
and of the Drude weight when approaching half-filling. RMFT and VMC calcu- 
lations also explain the large superconducting gap and the small quasiparticle 
weight in the underdoped Cuprates. Further, the modelling of charge ordered 
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states, impurity sites, or vortex cores qualitative agrees with experiments. 

While the quasiparticle weight Z vanishes in the half- filled limit, the nodal 
Fermi velocity vp stays finite. RMFT and VMC calculations explain this inter- 
esting experimental observation by the eflFect of the superexchange interaction 
on the dispersion. However, in the half-filled limit, such a behaviour {Z ^ 
and vp = const) immediately results in a divergence of the w- as well as of the 
k-dcpendence of the self-energy. The consequences of these divergences for any 
Fermi liquid description at finite doping are not fully understood. In a recent 
paper, Anderson [243] suggested that projected wave functions contain the es- 
sential physics to explain the non-Fermi liquid behaviour of the normal state 
in the cuprate superconductors, i.e. the region in the cuprate phase diagram 
above the pseudogap temperature scale. One reason why there has not been 
much progress on this issue is that we need a scheme to calculate the single 
particle Green's function directly in a Gutzwiller projected state. The standard 
technique of introducing a complete set of orthogonal excitations works as long 
as we only consider the contribution of the "projected quasiparticle (hole)" . 
However, as we discussed earlier, the total spectral weight of a photohole (say) 
is not exhausted by the projected quasihole excitation. The effect can be un- 
derstood most transparently as the non commutativity of a photohole state 
Cio-^bl^o) and projected excitations of the form PoQcrl^'o)- It follows that a 
photohole is a mixture of a projected single hole excitation and a multiparticle 
excitation which signifies the backflow of say,| spins accompanying a propagat- 
ing 1 spin hole. A consistent scheme to treat this effect has not been devised 
yet. 

Another important open question is the role played by phase fluctuations 
in Gutzwiller projected BCS wave functions. It was noted in the early papers 
of Anderson and collaborators that phase fluctuations are expected to play an 
important role as one approaches the Mott insulator in the phase diagram. The 
recent experiments of Ong and collaborator point to the existence of a vortex 
liquid phase above in the underdoped superconductors. A description of 
this phase within the Gutzwiller RVB theory has not yet been formulated. It 
should be pointed out that a large corpus of literature exists on fluctuating d 
- wave superconductors, but to our knowledge, no one has attempted to de- 
rive an effective "phase - only" model from a microscopic Hamiltonian for Mott 
Hubbard superconductors. Miiller-Hartmann and collaborators addressed this 
issue many years ago, using a high temperature expansion [244j . Using standard 
functional integral techniques, these authors performed a Ginzburg Landau ex- 
pansion for the free energy functional of an RVB state |130] . They flnd that 
singlet pairing sets in at temperature scales higher than the (mean field) tran- 
sition temperature T^. As hole concentration goes to zero, the local U{1) gauge 
symmetry in the theory leads to phase fluctuations that destroy off diagonal 
long range order. However, their calculations show that an extended s - wave 
state is favoured, which result is inconsistent with the Gutzwiller RVB theory 
described in this article. It would be very interesting to revisit this problem 
and attempt to derive a Ginzburg Landau expansion of the d - wave Gutzwiller 
RVB state. Such a step is necessary to extend the Gutzwiller - RVB framework 
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to the description of phase degrees of freedom, and the effect of the latter in 
destroying superconducting correlations, both as a function of temperature and 
doping. 

A related issue is the understanding of the pseudogap state within the RVB 
theory. The view we advocated was that local singlet pairing exists at temper- 
ature scales T < T*. Much support for this idea comes from the experimental 
observation that the BCS ratio, A/(2fcBT*), is constant and in agreement with 
mean field theory for all doping levels, when we use the onset temperature of 
the pseudogap T* instead of Tc- While this is certainly suggestive, there is no 
direct way of proving this within the theory, again because we do not yet have a 
method to describe finite temperature effects within the Gutzwiller framework. 
Extending the Gutzwiller-RVB theory to the description of finite temperature 
phases is an important step that needs to be taken to complete our understand- 
ing of the pseudogap state. In this context, we note that a a finite temperature 
scheme for Gutzwiller projected Fermi liquids was developed by Seller et al. to 
study '^He |121| . Whether a similar scheme can be developed to study projected 
d - wave superconductors at finite temperatures remains to be seen. A related 
issue is the investigation of finite frequency excitations in the Gutzwiller - RVB 
scheme. In particular, it will be very useful to study the collective excitations 
of the Gutzwiller superconductor along the lines of Anderson's original work on 
equations of motion for collective modes in a BCS superconductor |245j . 
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